1. Description

This R Markdown document describes the analyses performed for the manuscript entitled “Environmental pollution correlates with parasite infection across a riverine landscape” by Io S. Deflem, Seppe Marchand, Federico C.F. Calboli, Joost A.M. Raeymaekers, Filip A.M. Volckaert and Pascal I. Hablützel.

The analyses were run in R 4.2.2

2. Study area and sampling

Up to thirty 0+ three-spined sticklebacks were sampled at 37 locations in the Dijle and Demer basins in Belgium during autumn 2016 under a permit of the Flemish Agency Nature and Forest. Both basins together cover a continuous surface area of 3,624 km² with the furthest two sampling sites being located 117 km apart (distance measured along rivers). All locations included small and relatively slow flowing streams (drop off from highest to lowest point is 18 m) covering a wide range of ecological, hydromorphological, and physico-chemical characteristics. Fish were caught using a dip net.

3. Setting up working environment

# Empty environment
rm(list=ls())

# Set working directory to location where script is stored
setwd(dirname(rstudioapi::getActiveDocumentContext()$path)) # requires installation of package "rstudioapi"

# Loading required libraries
require(BAS)
require(boral)
require(car)
require(corrplot)
require(ggplot2)
require(gplots)
require(vegan)

4 Loading and preparing host and parasite data

Fish were euthanized with a lethal dose of MS222 on the day of capture, following directions of the KU Leuven Animal Ethics Commission, and stored at -20 °C. Fish were kept in separate containers per site at all times. Lab based parasite screening of thawed fish involved placing individual fish in 5 or 10 ml cryo-tubes with 1 or 2 ml of distilled water. Following a vigorous shake of 10 s, the liquid was poured into a Petri dish and ectoparasites were identified and counted using a stereomicroscope. Fish were rinsed and checked again for the presence of ectoparasites on skin and fins. The intestines were examined for endoparasites. Before dissection, fish weight (± 0.1 mg) and standard length (± 1 mm) were recorded. To quantify body condition, we calculated the scaled mass index (SMI; Maceda-Veiga et al., 2014; Peig & Green, 2009). Sex was determined during dissection by inspection of gonad development. A total of 668 fish were dissected, which amounts to approximately 20 fish per location, with the exception of seven locations where only 10 fish were screened for the presence of macroparasites. Ecto- and endoparasites were morphologically identified to species level whenever possible.

# Parasite data
data <- read.csv("data_2016_2303.csv", sep=';')
data$site <- as.factor(data$site)

# Calculate parasite parameters
#names(data)
#parasite data is overdispersed (mostly so for Trichodina), if using average abundance data, species matrix needs to be transformed
datao <- na.omit(data[,c(1,22:24,26:32)]) #remove fish without parasite counts

ddata <- dispweight(datao[,-1]) #correct for overdispersion of the parasite count data
avab <- aggregate(ddata, by = list(datao[,1]), function(x){mean(x, na.rm =T)})
prev = aggregate(data[,c(22:24,26:32)], by = list(data[,1]), function(x){sum(x >0, na.rm = T)/length(x)})
medin = aggregate(data[,c(22:24,26:32)], by = list(data[,1]), function(x){median(x[x >0], na.rm = T)}) 
pa = aggregate(data[,c(22:24,26:32)], by = list(data[,1]), function(x){ifelse(mean(x, na.rm =T)>0,1,0)}) 

avab[is.na(avab)] <- 0
prev[is.na(prev)] <- 0
medin[is.na(medin)] <- 0

# Host condition data
avcondition <- aggregate(data$SMI, by = list(data[,1]), function(x){mean(x, na.rm =T)})[,2]
avlength <- aggregate(data$length, by = list(data[,1]), function(x){mean(x, na.rm =T)})[,2]

# Parasite index
sgyr <- 1:nrow(data)
stri <- 1:nrow(data)
sglu <- 1:nrow(data)
scon <- 1:nrow(data)
scysl <- 1:nrow(data)
spro <- 1:nrow(data)
saca <- 1:nrow(data)
scam <- 1:nrow(data)
sang <- 1:nrow(data)

for(j in 1:nrow(data)){
  sgyr[j] <- data$Gyr[j]/sd(data$Gyr, na.rm=T)
  stri[j] <- data$Tri[j]/sd(data$Tri, na.rm=T)
  sglu[j] <- data$Glu[j]/sd(data$Glu, na.rm=T)
  scon[j] <- data$Con[j]/sd(data$Con, na.rm=T)
  scysl[j] <- data$CysL[j]/sd(data$CysL, na.rm=T)
  spro[j] <- data$Pro[j]/sd(data$Pro, na.rm=T)
  saca[j] <- data$Aca[j]/sd(data$Aca, na.rm=T)
  scam[j] <- data$Cam[j]/sd(data$Cam, na.rm=T)
  sang[j] <- data$Ang[j]/sd(data$Ang, na.rm=T)
}

PI <- 1:nrow(data)
for(j in 1:nrow(data)){
  PI[j] <- 10/max(sgyr, na.rm=T)*sgyr[j] + 10/max(stri, na.rm=T)*stri[j] + 10/max(sglu, na.rm=T)*sglu[j] +
    10/max(scon, na.rm=T)*scon[j] + 10/max(scysl, na.rm=T)*scysl[j] + 10/max(spro, na.rm=T)*spro[j] +
    10/max(saca, na.rm=T)*saca[j] + 10/max(scam, na.rm=T)*scam[j] + 10/max(sang, na.rm=T)*sang[j]
}

PI_ecto <- 1:nrow(data)
for(j in 1:nrow(data)){
    PI_ecto[j] <- 10/max(sgyr, na.rm=T)*sgyr[j] + 10/max(stri, na.rm=T)*stri[j] + 10/max(sglu, na.rm=T)*sglu[j]
}

PI_endo <- 1:nrow(data)
for(j in 1:nrow(data)){
  PI_endo[j] <-     10/max(scon, na.rm=T)*scon[j] + 10/max(scysl, na.rm=T)*scysl[j] + 10/max(spro, na.rm=T)*spro[j] +
    10/max(saca, na.rm=T)*saca[j] + 10/max(scam, na.rm=T)*scam[j] + 10/max(sang, na.rm=T)*sang[j]
}

avPI <- aggregate(PI, by = list(data[,1]), function(x){mean(x, na.rm =T)})[,2]
avPI_ecto <- aggregate(PI_ecto, by = list(data[,1]), function(x){mean(x, na.rm =T)})[,2]
avPI_endo <- aggregate(PI_endo, by = list(data[,1]), function(x){mean(x, na.rm =T)})[,2]

5 Loading and preparing environmental and spatial data

Physico-chemical data was provided by the Flemish Environmental Agency (VMM). Each fish sampling site was chosen at or near an environmental monitoring site of VMM. Water parameters include water temperature, pH, conductivity, dissolved oxygen (O2), saturation with dissolved oxygen, and Biochemical and Chemical Oxygen Demand (BOD and COD). Nutrient related water parameters include measurements of nitrate (NO3-), nitrite (NO2-), Kjeldahl nitrogen (KjN), total nitrogen (Nt), ammonium (NH4+), and total phosphorus (Pt). Following removal of strong collinear variables (significant correlation with P < 0.05 and Pearson correlation coefficient > |0.6|; Dormann et al., 2013), six environmental physico-chemical variables were retained (temperature, conductivity, COD, saturation with dissolved oxygen, ammonium, and total nitrogen), representing different aspects of water quality. For each water parameter, the average value of the year before sampling was calculated based on monthly monitoring data. Additionally, two hydromorphological variables were included: Tthe presence of a pool-riffle pattern and meanders were noted during field sampling and these parameters were included as binary variables (presence/absence) for a representation of abiotic habitat structure. Spatial (waterway) distances were calculated using the Network Analyst toolbox in ArcGIS. Upstream distance was defined as the maximal upstream distance from the sampling location. Network peripherality was calculated as the average waterway distance of a sampling location to all other locations. Hence, a total of eight environmental and two spatial variables were included in the statistical analysis.

# Environmental data (VMM)
environment <- read.csv("Environment_update.csv", sep=';')
spavar <- read.csv("space2.csv", sep=';') #spatial variables: network centrality and upstream distance
plot(spavar$netcen); plot(density(spavar$netcen))

plot(spavar$updist); plot(density(spavar$updist))

# Environmental data (from field observations)
field_data <- read.csv("field_data.csv", sep=',')
environment2 <- cbind(environment[,c(1,49,52:53,55,57,60,63)], field_data[-c(8,10,25,27,37),33:34], spavar[,c(2,3)])
environment2$pool_riffle <- as.factor(environment2$pool_riffle)
environment2$meander <- as.factor(environment2$meander)

netcen <- spavar$netcen
updist <- spavar$updist

We used univariate generalized linear models to investigate how landscape-level effects modify infection patterns of individual parasite taxa, host size and condition. We kept the statistical models linear (as opposed to polynomial) and only considered main effects (i.e. no interaction terms) because we had no prior information from this study system that more complex models were necessary and because the study design with (only) 37 sampling sites was not intended for non-linear models. Ten explanatory variables (temperature, conductivity, COD, saturation with dissolved oxygen, ammonium, total nitrogen, the presence of pool-riffle patterns and meanders, upstream distance, and network peripherality) were included.

6. Univeriate analysis using Bayesian Model Averaging

Univariate analyses - We used generalized linear models in a BMA approach to understand how infection with individual parasite taxa relate to host characteristics (length and condition), environmental conditions as well as spatial distance among sampling sites. Parasite infection was calculated in three ways at the host population level: average abundance (mean parasites per host), prevalence (percentage of infected hosts) and median infection intensity (median number of parasites in infected hosts). We calculated the individual parasitation index (IPI) following Kalbe et al. (2002) as a measurement for total parasite abundance and species richness for each individual fish. This index was calculated for all parasite species combined, and for ecto- and endoparasite species separately. For these models, we assumed a normal error distribution (which appeared to be justified, see Supplementary Figures S1-S2) and applied a Jeffrey-Zellner-Siow prior. Model assumptions (homoscedasticity of the variances and normal distribution of the errors) were assessed using the generic model plot function in R and did onlynot clearly deviate in any of the models for rare parasites. We followed a normal distribution, and not Poisson or negative binomial, for the parasite data for the common species (Trichodina sp. and Gyrodactylus spp.) and the individual parasitation index as the parameters used are deviates from count data. Rare parasites (Glugea, Contracaecum, Anguillicoloides, and unidentified cysts) were excluded from the univariate analysis because there was not enough data to obtain a good fit of the models.For rare parasites (Contracaecum sp. and Anguillicoloides crassus), we used population-level presence-absence data assuming a binomial error distribution and a uniformly distributed BIC prior. Due to low prevalences, the other parasites were not included in the species-specific models. Explanatory variables were considered important when they had a posterior inclusion probability (PIP) of 0.5. To account for overdispersion in the parasite counts, we transformed the data by downweighting overdispersed taxa following Clarke et al. (2006) using the dispweight function in the R package vegan v2.5.6 (Oksanen et al., 2013).

6.0.1 Figure 2

# Make a matrix for PIP (Posterior Inclusion Probability)
PIP <- matrix(nrow=12, ncol=14)
rownames(PIP) <- c("Host condition", "Host length", "Temperature", "Oxygen saturation", "Conductivity", "COD", "Ammonium", "Total nitrogen", "Pool riffle pattern", "Meander", "Network peripherality", "Upstream distance")
colnames(PIP) <- c("Host condition", "Host size", "Gyrodactylus abundance", "Gyrodactylus prevalence", "Gyrodactylus infection intensity", "Trichodina abundance", "Trichodina prevalence", "Trichodina infection intensity", "Glugea", "Contracaecum", "Aguillicola",
                   "PI", "PI ecto", "PI endo")  
#Condition
bas.model <- bas.lm(avcondition ~ T_av + O2_sat_av + Con_av + COD_av 
                     + NH4._av + Nt_av + pool_riffle + meander + netcen + 
                       updist, data=environment2, prior="JZS")
coef.model <- coef(bas.model)
pip <- summary(bas.model)
PIP[c(3:12),1] <- pip[2:11,1]*sign(coef.model$postmean[2:11])

#Length
bas.model <- bas.lm(avlength ~ T_av + O2_sat_av + Con_av + COD_av 
                     + NH4._av + Nt_av + pool_riffle + meander + netcen + 
                       updist, data=environment2, prior="JZS")
coef.model <- coef(bas.model)
pip <- summary(bas.model)
PIP[c(3:12),2] <- pip[2:11,1]*sign(coef.model$postmean[2:11])

#Gyrodactylus abundance
bas.model <- bas.lm(avab$Gyr ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av 
                     + NH4._av + Nt_av + pool_riffle + meander + netcen + 
                       updist, data=environment2, prior="JZS")
coef.model <- coef(bas.model)
pip <- summary(bas.model)
PIP[c(1:12),3] <- pip[2:13,1]*sign(coef.model$postmean[2:13])

#Gyrodactylus prevalence
bas.model <- bas.lm(prev$Gyr ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av 
                     + NH4._av + Nt_av + pool_riffle + meander + netcen + 
                       updist, data=environment2, prior="JZS")
coef.model <- coef(bas.model)
pip <- summary(bas.model)
PIP[c(1:12),4] <- pip[2:13,1]*sign(coef.model$postmean[2:13])

#Gyrodactylus infection intensity
bas.model <- bas.lm(medin$Gyr ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av 
                     + NH4._av + Nt_av + pool_riffle + meander + netcen + 
                       updist, data=environment2, prior="JZS")
coef.model <- coef(bas.model)
pip <- summary(bas.model)
PIP[c(1:12),5] <- pip[2:13,1]*sign(coef.model$postmean[2:13])

#Trichodina abundance
bas.model <- bas.lm(avab$Tri ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av 
                     + NH4._av + Nt_av + pool_riffle + meander + netcen + 
                       updist, data=environment2, prior="JZS")
coef.model <- coef(bas.model)
pip <- summary(bas.model)
PIP[c(1:12),6] <- pip[2:13,1]*sign(coef.model$postmean[2:13])

#Trichodina prevalence
bas.model <- bas.lm(prev$Tri ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av 
                     + NH4._av + Nt_av + pool_riffle + meander + netcen + 
                       updist, data=environment2, prior="JZS")
coef.model <- coef(bas.model)
pip <- summary(bas.model)
PIP[c(1:12),7] <- pip[2:13,1]*sign(coef.model$postmean[2:13])

#Trichodina infection intensity
bas.model <- bas.lm(medin$Tri ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av 
                     + NH4._av + Nt_av + pool_riffle + meander + netcen + 
                       updist, data=environment2, prior="JZS")
coef.model <- coef(bas.model)
pip <- summary(bas.model)
PIP[c(1:12),8] <- pip[2:13,1]*sign(coef.model$postmean[2:13])

#Glugea
bas.model <- bas.glm(pa$Glu ~  avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av 
                     + NH4._av + Nt_av + pool_riffle + meander + netcen + 
                       updist, data=environment2, betaprior=g.prior(100), family=binomial)
coef.model <- coef(bas.model)
pip <- summary(bas.model)
PIP[c(1:12),9] <- pip[2:13,1]*sign(coef.model$postmean[2:13])

#Contracaecum
bas.model <- bas.glm(pa$Con ~  avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av 
                     + NH4._av + Nt_av + pool_riffle + meander + netcen + 
                       updist, data=environment2, betaprior=g.prior(100), family=binomial)
coef.model <- coef(bas.model)
pip <- summary(bas.model)
PIP[c(1:12),10] <- pip[2:13,1]*sign(coef.model$postmean[2:13])

#Anguillicola
bas.model <- bas.glm(pa$Ang ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av 
                     + NH4._av + Nt_av + pool_riffle + meander + netcen + 
                       updist, data=environment2, betaprior=g.prior(100), family=binomial)
coef.model <- coef(bas.model)
pip <- summary(bas.model)
PIP[c(1:12),11] <- pip[2:13,1]*sign(coef.model$postmean[2:13])

#PI
bas.model <- bas.lm(avPI ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av 
                     + NH4._av + Nt_av + pool_riffle + meander + netcen + 
                       updist, data=environment2, prior="JZS")
coef.model <- coef(bas.model)
pip <- summary(bas.model)
PIP[c(1:12),12] <- pip[2:13,1]*sign(coef.model$postmean[2:13])

#PI ecto
bas.model <- bas.lm(avPI_ecto ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av 
                     + NH4._av + Nt_av + pool_riffle + meander + netcen + 
                       updist, data=environment2, prior="JZS")
coef.model <- coef(bas.model)
pip <- summary(bas.model)
PIP[c(1:12),13] <- pip[2:13,1]*sign(coef.model$postmean[2:13])

#PI endo
bas.model <- bas.lm(avPI_endo ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av 
                     + NH4._av + Nt_av + pool_riffle + meander + netcen + 
                       updist, data=environment2, prior="JZS")
coef.model <- coef(bas.model)
pip <- summary(bas.model)
PIP[c(1:12),14] <- pip[2:13,1]*sign(coef.model$postmean[2:13])
x = round(PIP, digits=2)
x[abs(PIP)<0.5] <- ""
x[abs(PIP)>0.5] <- "+"
heatmap.2(PIP[,-c(9,10,11)],
          cellnote = x[,-c(9,10,11)],
          #main = "Correlation",
          notecex=1,
          notecol="white",
          density.info="none",
          trace="none",
          margins =c(10,8),
          col=redblue(256),
          dendrogram="both",
          cexRow = 0.7,
          cexCol = 0.7,
          key.title = "PIP",
          lhei = c(1,3),
          lwid = c(0.5, 0.5),
          srtCol = 45,
          labCol = c("Host condition", "Host length", expression(paste(italic("Gyrodactylus"),  " spp. - abundance")), expression(paste(italic("Gyrodactylus"), " spp. - prevalence")), expression(paste(italic("Gyrodactylus"), " spp. - infection intensity")), expression(paste(italic("Trichodina"), " sp. - abundance")), expression(paste(italic("Trichodina"), " sp. - prevalence")), expression(paste(italic("Trichodina"), " sp. - infection intensity")), expression("I"[PI]*" all parasites"), expression("I"[PI]*" ectoparasites"), expression("I"[PI]*" endoparasites")) ,
          #Colv="NA"
          ) 

pdf("Figure2.pdf", width = 7.29, height = 4.5)
heatmap.2(PIP[,-c(9,10,11)],
          cellnote = x[,-c(9,10,11)],
          #main = "Correlation",
          notecex=1,
          notecol="white",
          density.info="none",
          trace="none",
          margins =c(10,8),
          col=redblue(256),
          dendrogram="both",
          cexRow = 0.7,
          cexCol = 0.7,
          key.title = "PIP",
          lhei = c(1,3),
          lwid = c(0.5, 0.5),
          srtCol = 45,
          labCol = c("Host condition", "Host length", expression(paste(italic("Gyrodactylus"),  " spp. - abundance")), expression(paste(italic("Gyrodactylus"), " spp. - prevalence")), expression(paste(italic("Gyrodactylus"), " spp. - infection intensity")), expression(paste(italic("Trichodina"), " sp. - abundance")), expression(paste(italic("Trichodina"), " sp. - prevalence")), expression(paste(italic("Trichodina"), " sp. - infection intensity")), expression("I"[PI]*" all parasites"), expression("I"[PI]*" ectoparasites"), expression("I"[PI]*" endoparasites")) ,
          #Colv="NA"
          ) 

dev.off()
## png 
##   2

6.1 Variation in host condition

bas.model <- bas.lm(avcondition ~  T_av + O2_sat_av + Con_av + COD_av + NH4._av + Nt_av + 
                      pool_riffle + meander + netcen + updist, 
                    data=environment2, prior="JZS")
plot(bas.model)

summary(bas.model)
##              P(B != 0 | Y) model 1    model 2    model 3    model 4    model 5
## Intercept       1.00000000  1.0000  1.0000000  1.0000000  1.0000000  1.0000000
## T_av            0.12386945  0.0000  1.0000000  0.0000000  0.0000000  0.0000000
## O2_sat_av       0.05073157  0.0000  0.0000000  0.0000000  0.0000000  0.0000000
## Con_av          0.03919858  0.0000  0.0000000  0.0000000  0.0000000  0.0000000
## COD_av          0.05197281  0.0000  0.0000000  0.0000000  0.0000000  1.0000000
## NH4._av         0.04846220  0.0000  0.0000000  0.0000000  0.0000000  0.0000000
## Nt_av           0.05248386  0.0000  0.0000000  0.0000000  0.0000000  0.0000000
## pool_riffle1    0.06729856  0.0000  0.0000000  0.0000000  1.0000000  0.0000000
## meander1        0.06527705  0.0000  0.0000000  0.0000000  0.0000000  0.0000000
## netcen          0.06898689  0.0000  0.0000000  1.0000000  0.0000000  0.0000000
## updist          0.04496936  0.0000  0.0000000  0.0000000  0.0000000  0.0000000
## BF                      NA  1.0000  0.6511099  0.3462387  0.2957714  0.2212089
## PostProbs               NA  0.6912  0.0450000  0.0239000  0.0204000  0.0153000
## R2                      NA  0.0000  0.0906000  0.0565000  0.0478000  0.0315000
## dim                     NA  1.0000  2.0000000  2.0000000  2.0000000  2.0000000
## logmarg                 NA  0.0000 -0.4290769 -1.0606268 -1.2181684 -1.5086476
image(bas.model, rotate=F)

coef.model <- coef(bas.model)
#abs(coef.model$postmean)-2*coef.model$postsd > 0
plot(confint(coef.model, parm = 2:11))
## Warning in arrows(x[not.deg], ci[not.deg, 1], x[not.deg], ci[not.deg, 2], :
## zero-length arrow is of indeterminate angle and so skipped

## NULL
confint <- confint(coef.model, parm = 2:11)
write.table(confint, 'condition.txt', sep="\t")
pip <- summary(bas.model)
PIP[c(3:12),1] <- pip[2:11,1]*sign(coef.model$postmean[2:11])
#coef.model$postmean[2:11]

6.2 Variation in host length

bas.model <- bas.lm(avlength ~ T_av + O2_sat_av + Con_av + COD_av + NH4._av + Nt_av + 
                      pool_riffle + meander + netcen + updist, 
                    data=environment2, prior="JZS")
plot(bas.model)

summary(bas.model)
##              P(B != 0 | Y)   model 1    model 2   model 3  model 4   model 5
## Intercept        1.0000000 1.0000000 1.00000000 1.0000000 1.000000 1.0000000
## T_av             0.1914598 0.0000000 0.00000000 0.0000000 0.000000 0.0000000
## O2_sat_av        0.1407247 0.0000000 0.00000000 0.0000000 0.000000 0.0000000
## Con_av           0.4373524 0.0000000 0.00000000 1.0000000 1.000000 0.0000000
## COD_av           0.1445232 0.0000000 0.00000000 0.0000000 0.000000 0.0000000
## NH4._av          0.2026061 0.0000000 0.00000000 0.0000000 0.000000 0.0000000
## Nt_av            0.1889236 0.0000000 0.00000000 0.0000000 0.000000 0.0000000
## pool_riffle1     0.2259504 0.0000000 0.00000000 0.0000000 0.000000 1.0000000
## meander1         0.3217534 0.0000000 0.00000000 0.0000000 1.000000 0.0000000
## netcen           0.6121912 1.0000000 0.00000000 0.0000000 0.000000 1.0000000
## updist           0.1480063 0.0000000 0.00000000 0.0000000 0.000000 0.0000000
## BF                      NA 0.8278694 0.07277547 0.2894406 1.000000 0.6785609
## PostProbs               NA 0.1582000 0.13910000 0.0553000 0.042500 0.0288000
## R2                      NA 0.2306000 0.00000000 0.1819000 0.314300 0.2982000
## dim                     NA 2.0000000 1.00000000 2.0000000 3.000000 3.0000000
## logmarg                 NA 2.4314765 0.00000000 1.3805712 2.620376 2.2325953
image(bas.model, rotate=F)

coef.model <- coef(bas.model)
#abs(coef.model$postmean)-2*coef.model$postsd > 0
plot(confint(coef.model, parm = 2:11))
## Warning in arrows(x[not.deg], ci[not.deg, 1], x[not.deg], ci[not.deg, 2], :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(x[not.deg], ci[not.deg, 1], x[not.deg], ci[not.deg, 2], :
## zero-length arrow is of indeterminate angle and so skipped

## NULL
confint <- confint(coef.model, parm = 2:11)
write.table(confint, 'length.txt', sep="\t")
pip <- summary(bas.model)
PIP[c(3:12),1] <- pip[2:11,1]*sign(coef.model$postmean[2:11])
#coef.model$postmean[2:11]

# Prediction plot
newdata = as.data.frame(cbind(rep(mean(environment2$T_av), 37), 
                rep(mean(environment2$O2_sat_av), 37),
                rep(mean(environment2$Con_av), 37),
                rep(mean(environment2$COD_av), 37),
                rep(mean(environment2$NH4._av), 37),
                rep(mean(environment2$Nt_av), 37),
                rep(1, 37),
                rep(1, 37),
                rep(mean(netcen), 37),
                rep(mean(updist), 37)))
colnames(newdata) <- c("T_av", "O2_sat_av", "Con_av", "COD_av", "NH4._av", "Nt_av", "pool_riffle", "meander", "netcen", "updist")
newdata[,"pool_riffle"] <- as.factor(newdata[,"pool_riffle"]); newdata[,"meander"] <- as.factor(newdata[,"meander"])
newdata1 <- newdata; newdata1[,"netcen"] <- netcen
BMA_avlength_netcen <- predict(bas.model, newdata = newdata1, estimator = "BMA", se.fit=TRUE)

figure3j <- ggplot(environment2, aes(netcen, BMA_avlength_netcen$fit)) +
  theme_bw() +
  geom_line(color="red", size=1) +
  geom_ribbon(aes(ymin = (BMA_avlength_netcen$fit-BMA_avlength_netcen$se.bma.fit), ymax = (BMA_avlength_netcen$fit+BMA_avlength_netcen$se.bma.fit)), alpha = .1) +
  geom_point(data = environment2, aes(x=netcen, y=avlength)) +
  labs(x=expression("Network peripherality [m]"), y=expression("Average host length [mm]")) +
  theme(axis.title.x = element_text(size=12),
        axis.title.y = element_text(size=12)) +  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  labs(subtitle = "j)")
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
figure3j

6.3 Variation in Gyrodactylus infection

6.3.1 Mean abundance

bas.model <- bas.lm(avab$Gyr ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av 
                     + NH4._av + Nt_av + pool_riffle + meander + netcen + 
                       updist, data=environment2, prior="JZS")
yhat = fitted(bas.model, estimator = "BMA") #these are the fitted values under BMA
r = bas.model$Y - yhat #these are the model residuals
plot(bas.model)

summary(bas.model)
##              P(B != 0 | Y)  model 1    model 2   model 3   model 4   model 5
## Intercept        1.0000000 1.000000 1.00000000 1.0000000 1.0000000 1.0000000
## avlength         0.5500651 1.000000 0.00000000 0.0000000 0.0000000 1.0000000
## avcondition      0.3049274 0.000000 0.00000000 0.0000000 0.0000000 1.0000000
## T_av             0.1620852 0.000000 0.00000000 0.0000000 0.0000000 0.0000000
## O2_sat_av        0.2017699 0.000000 0.00000000 0.0000000 0.0000000 0.0000000
## Con_av           0.2621665 0.000000 0.00000000 0.0000000 0.0000000 0.0000000
## COD_av           0.8086300 1.000000 0.00000000 1.0000000 1.0000000 1.0000000
## NH4._av          0.2417809 0.000000 0.00000000 0.0000000 0.0000000 0.0000000
## Nt_av            0.8618345 1.000000 1.00000000 1.0000000 1.0000000 1.0000000
## pool_riffle1     0.1760484 0.000000 0.00000000 0.0000000 0.0000000 0.0000000
## meander1         0.1775842 0.000000 0.00000000 0.0000000 0.0000000 0.0000000
## netcen           0.3006370 0.000000 0.00000000 1.0000000 0.0000000 0.0000000
## updist           0.1648814 0.000000 0.00000000 0.0000000 0.0000000 0.0000000
## BF                      NA 1.000000 0.03902505 0.5452757 0.1537678 0.8271829
## PostProbs               NA 0.072100 0.05160000 0.0393000 0.0369000 0.0265000
## R2                      NA 0.524000 0.28790000 0.5059000 0.4100000 0.5631000
## dim                     NA 4.000000 2.00000000 4.0000000 3.0000000 5.0000000
## logmarg                 NA 6.997337 3.75378552 6.3908733 5.1250253 6.8076077
image(bas.model, rotate=F)

coef.model <- coef(bas.model)
abs(coef.model$postmean)-2*coef.model$postsd > 0
##  [1]  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [13] FALSE
confint(coef.model)
##                       2.5%        97.5%          beta
## Intercept     9.306741e-02 2.009914e-01  1.489134e-01
## avlength     -2.516824e-02 0.000000e+00 -8.122675e-03
## avcondition  -1.008942e+00 4.722253e-04 -1.500326e-01
## T_av         -1.933113e-02 1.936668e-02 -1.493132e-04
## O2_sat_av    -1.010761e-03 5.405876e-03  3.794809e-04
## Con_av       -7.801740e-05 5.075870e-04  5.726233e-05
## COD_av        0.000000e+00 1.295780e-02  6.652983e-03
## NH4._av      -6.043975e-02 7.670275e-02 -1.888173e-04
## Nt_av         0.000000e+00 5.253835e-02  2.930918e-02
## pool_riffle1 -3.248206e-02 1.064209e-01  3.966550e-03
## meander1     -9.521821e-02 4.524471e-02 -4.806259e-03
## netcen       -4.356805e-07 9.961756e-06  1.332021e-06
## updist       -1.428611e-06 1.623940e-06  4.618571e-08
## attr(,"Probability")
## [1] 0.95
## attr(,"class")
## [1] "confint.bas"
plot(confint(coef.model, parm = 2:11))

## NULL
confint <- confint(coef.model, parm = 2:11)
write.table(confint, 'GyroAA.txt', sep="\t")
pip <- summary(bas.model)
PIP[c(1:12),3] <- pip[2:13,1]*sign(coef.model$postmean[2:13])

6.3.1.1 Prediction plot for marginal effect of host condition on average Gyrodactylus infection

# Prediction plot
newdata = as.data.frame(cbind(rep(mean(avlength), 37),
                rep(mean(avcondition), 37),
                rep(mean(environment2$T_av), 37), 
                rep(mean(environment2$O2_sat_av), 37),
                rep(mean(environment2$Con_av), 37),
                rep(mean(environment2$COD_av), 37),
                rep(mean(environment2$NH4._av), 37),
                rep(mean(environment2$Nt_av), 37),
                rep(1, 37),
                rep(1, 37),
                rep(mean(netcen), 37),
                rep(mean(updist), 37)))
colnames(newdata) <- c("avlength", "avcondition", "T_av", "O2_sat_av", "Con_av", "COD_av", "NH4._av", "Nt_av", "pool_riffle", "meander", "netcen", "updist")
newdata[,"pool_riffle"] <- as.factor(newdata[,"pool_riffle"]); newdata[,"meander"] <- as.factor(newdata[,"meander"])
newdata1 <- newdata; newdata1[,"avcondition"] <- avcondition
BMA_Gyr_avcond <- predict(bas.model, newdata = newdata1, estimator = "BMA", se.fit=TRUE)
figure3k = ggplot(environment2, aes(avcondition, BMA_Gyr_avcond$fit)) +
  theme_bw() +
  geom_line(color="red", size=1) +
  geom_ribbon(aes(ymin = (BMA_Gyr_avcond$fit-BMA_Gyr_avcond$se.bma.fit), ymax = (BMA_Gyr_avcond$fit+BMA_Gyr_avcond$se.bma.fit)), alpha = .1) +
  geom_point(data = environment2, aes(x=avcondition, y=avab$Gyr)) +
  labs(x=expression("Average host condition"), y=expression(paste(italic("Gyrodactylus"), " spp. - abundance"))) +
  theme(axis.title.x = element_text(size=12),
        axis.title.y = element_text(size=12)) +  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  labs(subtitle = "k)")

figure3k

6.3.1.2 Prediction plot for marginal effect of COD on average Gyrodactylus infection

newdata1 <- newdata; newdata1[,"COD_av"] <- environment2$COD_av
BMA_Gyr_COD_av <- predict(bas.model, newdata = newdata1, estimator = "BMA", se.fit=TRUE)
figure3a = ggplot(environment2, aes(COD_av, BMA_Gyr_COD_av$fit)) +
  theme_bw() +
  geom_line(color="red", size=1) +
  geom_ribbon(aes(ymin = (BMA_Gyr_COD_av$fit-BMA_Gyr_COD_av$se.bma.fit), ymax = (BMA_Gyr_COD_av$fit+BMA_Gyr_COD_av$se.bma.fit)), alpha = .1) +
  geom_point(data = environment2, aes(x=COD_av, y=avab$Gyr)) +
  labs(x=expression("Chemical Oxygen Demand [mg/L]"), y=expression(paste(italic("Gyrodactylus"), " spp. - abundance"))) +
  theme(axis.title.x = element_text(size=12),
        axis.title.y = element_text(size=12)) +  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  labs(subtitle = "a)")
figure3a

6.3.2.3 Prediction plot for marginal effect of total nitrogen on average Gyrodactylus infection

newdata1 <- newdata; newdata1[,"Nt_av"] <- environment2$Nt_av
BMA_Gyr_Nt_av <- predict(bas.model, newdata = newdata1, estimator = "BMA", se.fit=TRUE)
figure3e = ggplot(environment2, aes(Nt_av, BMA_Gyr_Nt_av$fit)) +
  theme_bw() +
  geom_line(color="red", size=1) +
  geom_ribbon(aes(ymin = (BMA_Gyr_Nt_av$fit-BMA_Gyr_Nt_av$se.bma.fit), ymax = (BMA_Gyr_Nt_av$fit+BMA_Gyr_Nt_av$se.bma.fit)), alpha = .1) +
  geom_point(data = environment2, aes(x=Nt_av, y=avab$Gyr)) +
  labs(x=expression("Nt [mg/L]"), y=expression(paste(italic("Gyrodactylus"), " spp. - abundance"))) +
  theme(axis.title.x = element_text(size=12),
        axis.title.y = element_text(size=12)) +  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  labs(subtitle = "e)")
figure3e

6.3.2 Median infection intensity

bas.model <- bas.lm(medin$Gyr ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av 
                     + NH4._av + Nt_av + pool_riffle + meander + netcen + 
                       updist, data=environment2, prior="JZS")
yhat = fitted(bas.model, estimator = "BMA") #these are the fitted values under BMA
r = bas.model$Y - yhat #these are the model residuals

plot(bas.model)

summary(bas.model)
##              P(B != 0 | Y) model 1    model 2   model 3    model 4    model 5
## Intercept       1.00000000  1.0000  1.0000000  1.000000  1.0000000  1.0000000
## avlength        0.02065469  0.0000  0.0000000  0.000000  0.0000000  0.0000000
## avcondition     0.04286905  0.0000  1.0000000  0.000000  0.0000000  0.0000000
## T_av            0.02991678  0.0000  0.0000000  0.000000  0.0000000  0.0000000
## O2_sat_av       0.02400544  0.0000  0.0000000  0.000000  0.0000000  0.0000000
## Con_av          0.02143314  0.0000  0.0000000  0.000000  0.0000000  0.0000000
## COD_av          0.03655302  0.0000  0.0000000  0.000000  1.0000000  0.0000000
## NH4._av         0.03250611  0.0000  0.0000000  0.000000  0.0000000  1.0000000
## Nt_av           0.04088810  0.0000  0.0000000  1.000000  0.0000000  0.0000000
## pool_riffle1    0.02307275  0.0000  0.0000000  0.000000  0.0000000  0.0000000
## meander1        0.02038443  0.0000  0.0000000  0.000000  0.0000000  0.0000000
## netcen          0.02173649  0.0000  0.0000000  0.000000  0.0000000  0.0000000
## updist          0.02738794  0.0000  0.0000000  0.000000  0.0000000  0.0000000
## BF                      NA  1.0000  0.3176354  0.302852  0.2488803  0.2411373
## PostProbs               NA  0.7775  0.0206000  0.019600  0.0161000  0.0156000
## R2                      NA  0.0000  0.0517000  0.049100  0.0381000  0.0363000
## dim                     NA  1.0000  2.0000000  2.000000  2.0000000  2.0000000
## logmarg                 NA  0.0000 -1.1468512 -1.194511 -1.3907833 -1.4223890
image(bas.model, rotate=F)

coef.model <- coef(bas.model)
abs(coef.model$postmean)-2*coef.model$postsd > 0
##  [1]  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [13] FALSE
confint(coef.model)
##                  2.5%    97.5%          beta
## Intercept    1.607225 3.263557  2.405405e+00
## avlength     0.000000 0.000000 -5.938396e-04
## avcondition  0.000000 0.000000 -2.939823e-01
## T_av         0.000000 0.000000 -7.368692e-03
## O2_sat_av    0.000000 0.000000 -4.581964e-04
## Con_av       0.000000 0.000000 -1.594570e-05
## COD_av       0.000000 0.000000  1.609035e-03
## NH4._av      0.000000 0.000000  1.027436e-02
## Nt_av        0.000000 0.000000  7.570653e-03
## pool_riffle1 0.000000 0.000000  1.051497e-02
## meander1     0.000000 0.000000  2.511911e-03
## netcen       0.000000 0.000000  3.990331e-07
## updist       0.000000 0.000000  4.560731e-07
## attr(,"Probability")
## [1] 0.95
## attr(,"class")
## [1] "confint.bas"
plot(confint(coef.model, parm = 2:11))

## NULL
confint <- confint(coef.model, parm = 2:11)
write.table(confint, 'GyroAA.txt', sep="\t")
pip <- summary(bas.model)
PIP[c(1:12),3] <- pip[2:13,1]*sign(coef.model$postmean[2:13])

6.3.2 Prevalence

bas.model <- bas.lm(prev$Gyr ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av 
                     + NH4._av + Nt_av + pool_riffle + meander + netcen + 
                       updist, data=environment2, prior="JZS")
yhat = fitted(bas.model, estimator = "BMA") #these are the fitted values under BMA
r = bas.model$Y - yhat #these are the model residuals

plot(bas.model)

summary(bas.model)
##              P(B != 0 | Y)  model 1    model 2   model 3   model 4     model 5
## Intercept       1.00000000 1.000000 1.00000000 1.0000000 1.0000000  1.00000000
## avlength        0.16685753 0.000000 0.00000000 0.0000000 1.0000000  0.00000000
## avcondition     0.07728991 0.000000 0.00000000 0.0000000 0.0000000  0.00000000
## T_av            0.07242807 0.000000 0.00000000 0.0000000 0.0000000  0.00000000
## O2_sat_av       0.09667222 0.000000 0.00000000 0.0000000 0.0000000  0.00000000
## Con_av          0.58076129 1.000000 0.00000000 1.0000000 0.0000000  0.00000000
## COD_av          0.17268847 0.000000 0.00000000 1.0000000 0.0000000  0.00000000
## NH4._av         0.09134534 0.000000 0.00000000 0.0000000 0.0000000  0.00000000
## Nt_av           0.11056039 0.000000 0.00000000 0.0000000 0.0000000  0.00000000
## pool_riffle1    0.13922799 0.000000 0.00000000 0.0000000 0.0000000  0.00000000
## meander1        0.09905733 0.000000 0.00000000 0.0000000 0.0000000  0.00000000
## netcen          0.11864979 0.000000 0.00000000 0.0000000 0.0000000  1.00000000
## updist          0.06837623 0.000000 0.00000000 0.0000000 0.0000000  0.00000000
## BF                      NA 1.000000 0.08283414 0.8857497 0.1246735  0.07409002
## PostProbs               NA 0.228800 0.22750000 0.0369000 0.0285000  0.01700000
## R2                      NA 0.233300 0.00000000 0.3039000 0.1341000  0.10730000
## dim                     NA 2.000000 1.00000000 3.0000000 2.0000000  2.00000000
## logmarg                 NA 2.490915 0.00000000 2.3695941 0.4088580 -0.11155941
image(bas.model, rotate=F)

coef.model <- coef(bas.model)
abs(coef.model$postmean)-2*coef.model$postsd > 0
##  [1]  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [13] FALSE
confint(coef.model)
##                       2.5%        97.5%          beta
## Intercept     2.547614e-01 3.964116e-01  3.224082e-01
## avlength     -2.300336e-02 0.000000e+00 -2.431148e-03
## avcondition  -3.425493e-01 9.865432e-02 -1.986363e-02
## T_av         -1.891814e-02 2.173153e-04 -7.123130e-04
## O2_sat_av    -4.068781e-03 1.262431e-04 -2.662986e-04
## Con_av        0.000000e+00 8.064140e-04  3.112530e-04
## COD_av        0.000000e+00 8.036852e-03  9.376253e-04
## NH4._av      -3.593614e-03 4.398906e-02  1.064409e-03
## Nt_av        -1.789721e-04 2.438631e-02  1.612509e-03
## pool_riffle1  0.000000e+00 1.577918e-01  1.465705e-02
## meander1     -1.043718e-01 4.991691e-04 -6.966751e-03
## netcen       -2.974468e-09 8.169806e-06  5.765326e-07
## updist       -8.239812e-07 1.990176e-07 -6.922933e-09
## attr(,"Probability")
## [1] 0.95
## attr(,"class")
## [1] "confint.bas"
plot(confint(coef.model, parm = 2:11))

## NULL
confint <- confint(coef.model, parm = 2:11)
write.table(confint, 'GyroAA.txt', sep="\t")
pip <- summary(bas.model)
PIP[c(1:12),3] <- pip[2:13,1]*sign(coef.model$postmean[2:13])

6.3.2.3 Prediction plot for marginal effect of conductivity on average Gyrodactylus infection

newdata1 <- newdata; newdata1[,"Con_av"] <- environment2$Con_av
BMA_Gyr_Con_av <- predict(bas.model, newdata = newdata1, estimator = "BMA", se.fit=TRUE)
figure3i = ggplot(environment2, aes(Con_av, BMA_Gyr_Con_av$fit)) +
  theme_bw() +
  geom_line(color="red", size=1) +
  geom_ribbon(aes(ymin = (BMA_Gyr_Con_av$fit-BMA_Gyr_Con_av$se.bma.fit), ymax = (BMA_Gyr_Con_av$fit+BMA_Gyr_Con_av$se.bma.fit)), alpha = .1) +
  geom_point(data = environment2, aes(x=Con_av, y=prev$Gyr)) +
  labs(x=expression(paste("Conductivity [", mu, "S/cm]")), y=expression(paste(italic("Gyrodactylus"), " spp. - prevalence"))) +
  theme(axis.title.x = element_text(size=12),
        axis.title.y = element_text(size=12)) +  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  labs(subtitle = "i)")
figure3i

6.4 Variation in Trichodina infection

6.4.1 Mean abundance

bas.model <- bas.lm(avab$Tri ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av 
                     + NH4._av + Nt_av + pool_riffle + meander + netcen + 
                       updist, data=environment2, prior="JZS")
yhat = fitted(bas.model, estimator = "BMA") #these are the fitted values under BMA
r = bas.model$Y - yhat #these are the model residuals
plot(bas.model)

summary(bas.model)
##              P(B != 0 | Y)   model 1  model 2  model 3   model 4   model 5
## Intercept       1.00000000 1.0000000 1.000000 1.000000 1.0000000 1.0000000
## avlength        0.12444042 0.0000000 0.000000 0.000000 0.0000000 0.0000000
## avcondition     0.09559152 0.0000000 0.000000 0.000000 0.0000000 0.0000000
## T_av            0.08500185 0.0000000 0.000000 0.000000 0.0000000 0.0000000
## O2_sat_av       0.09751169 0.0000000 0.000000 0.000000 0.0000000 0.0000000
## Con_av          0.42180029 0.0000000 1.000000 0.000000 0.0000000 0.0000000
## COD_av          0.46724539 0.0000000 1.000000 0.000000 0.0000000 1.0000000
## NH4._av         0.21465122 0.0000000 0.000000 0.000000 1.0000000 0.0000000
## Nt_av           0.34576751 0.0000000 0.000000 1.000000 0.0000000 1.0000000
## pool_riffle1    0.13743699 0.0000000 0.000000 0.000000 0.0000000 0.0000000
## meander1        0.08968689 0.0000000 0.000000 0.000000 0.0000000 0.0000000
## netcen          0.10313787 0.0000000 0.000000 0.000000 0.0000000 0.0000000
## updist          0.08604777 0.0000000 0.000000 0.000000 0.0000000 0.0000000
## BF                      NA 0.0239404 1.000000 0.170643 0.1105942 0.2711811
## PostProbs               NA 0.1650000 0.104400 0.098000 0.0635000 0.0283000
## R2                      NA 0.0000000 0.358500 0.209300 0.1890000 0.3063000
## dim                     NA 1.0000000 3.000000 2.000000 2.0000000 3.0000000
## logmarg                 NA 0.0000000 3.732188 1.964006 1.5303004 2.4272192
image(bas.model, rotate=F)

coef.model <- coef(bas.model)
abs(coef.model$postmean)-2*coef.model$postsd > 0
##  [1]  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [13] FALSE
confint(coef.model)
##                       2.5%        97.5%          beta
## Intercept     1.232338e-01 2.714579e-01  2.002413e-01
## avlength     -1.705124e-02 3.780264e-04 -1.114390e-03
## avcondition  -5.882325e-01 5.136492e-02 -2.394139e-02
## T_av         -1.537747e-02 1.020332e-02 -1.479738e-04
## O2_sat_av    -1.352662e-03 2.779387e-03  1.132458e-04
## Con_av        0.000000e+00 7.626598e-04  2.109076e-04
## COD_av       -1.738195e-05 1.373697e-02  4.203136e-03
## NH4._av      -8.071803e-04 1.061136e-01  1.096330e-02
## Nt_av         0.000000e+00 4.926715e-02  1.044647e-02
## pool_riffle1  0.000000e+00 1.494898e-01  1.212839e-02
## meander1     -7.138422e-02 2.663792e-02 -1.678368e-03
## netcen       -6.434350e-06 4.469308e-07 -2.846749e-07
## updist       -1.450572e-06 7.187041e-07 -1.915234e-08
## attr(,"Probability")
## [1] 0.95
## attr(,"class")
## [1] "confint.bas"
plot(confint(coef.model, parm = 2:11))

## NULL
confint <- confint(coef.model, parm = 2:11)
write.table(confint, 'GyroAA.txt', sep="\t")
pip <- summary(bas.model)
PIP[c(1:12),3] <- pip[2:13,1]*sign(coef.model$postmean[2:13])

6.4.1 Median infection intensity

bas.model <- bas.lm(medin$Tri ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av 
                     + NH4._av + Nt_av + pool_riffle + meander + netcen + 
                       updist, data=environment2, prior="JZS")
yhat = fitted(bas.model, estimator = "BMA") #these are the fitted values under BMA
r = bas.model$Y - yhat #these are the model residuals
plot(bas.model)

summary(bas.model)
##              P(B != 0 | Y)  model 1    model 2   model 3   model 4   model 5
## Intercept        1.0000000 1.000000 1.00000000 1.0000000 1.0000000 1.0000000
## avlength         0.2881774 0.000000 0.00000000 0.0000000 0.0000000 0.0000000
## avcondition      0.2224163 0.000000 0.00000000 0.0000000 1.0000000 0.0000000
## T_av             0.1574615 0.000000 0.00000000 0.0000000 0.0000000 0.0000000
## O2_sat_av        0.1621652 0.000000 0.00000000 0.0000000 0.0000000 0.0000000
## Con_av           0.7912882 1.000000 0.00000000 1.0000000 1.0000000 1.0000000
## COD_av           0.7219333 1.000000 0.00000000 0.0000000 1.0000000 1.0000000
## NH4._av          0.3248706 0.000000 1.00000000 1.0000000 0.0000000 0.0000000
## Nt_av            0.3218716 0.000000 0.00000000 0.0000000 0.0000000 0.0000000
## pool_riffle1     0.2081402 0.000000 0.00000000 0.0000000 0.0000000 0.0000000
## meander1         0.2655250 0.000000 0.00000000 0.0000000 0.0000000 1.0000000
## netcen           0.2107764 0.000000 0.00000000 0.0000000 0.0000000 0.0000000
## updist           0.1536937 0.000000 0.00000000 0.0000000 0.0000000 0.0000000
## BF                      NA 1.000000 0.04957628 0.2279222 0.5092549 0.4837322
## PostProbs               NA 0.148800 0.04060000 0.0339000 0.0227000 0.0216000
## R2                      NA 0.449900 0.26810000 0.3987000 0.4816000 0.4799000
## dim                     NA 3.000000 2.00000000 3.0000000 4.0000000 4.0000000
## logmarg                 NA 6.288599 3.28435583 4.8098478 5.6137920 5.5623749
image(bas.model, rotate=F)

coef.model <- coef(bas.model)
abs(coef.model$postmean)-2*coef.model$postsd > 0
##  [1]  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [13] FALSE
confint(coef.model)
##                       2.5%        97.5%          beta
## Intercept     1.675518e+01 3.441505e+01  2.547297e+01
## avlength     -3.471106e+00 1.978629e-01 -5.158998e-01
## avcondition  -1.452798e+02 8.422586e-02 -1.375058e+01
## T_av         -1.669965e+00 4.431680e+00  2.065105e-01
## O2_sat_av    -4.995746e-01 4.824268e-01 -1.498026e-02
## Con_av        0.000000e+00 1.283962e-01  6.585727e-02
## COD_av        0.000000e+00 2.001430e+00  9.195418e-01
## NH4._av      -2.433043e+00 1.595274e+01  2.222211e+00
## Nt_av        -3.006566e-01 6.163537e+00  1.012166e+00
## pool_riffle1 -2.712539e+00 2.385458e+01  2.247962e+00
## meander1     -2.677652e+01 2.623704e-01 -3.523477e+00
## netcen       -1.459503e-03 1.344122e-04 -1.366046e-04
## updist       -2.059586e-04 2.662744e-04  5.653370e-06
## attr(,"Probability")
## [1] 0.95
## attr(,"class")
## [1] "confint.bas"
plot(confint(coef.model, parm = 2:11))

## NULL
confint <- confint(coef.model, parm = 2:11)
write.table(confint, 'GyroAA.txt', sep="\t")
pip <- summary(bas.model)
PIP[c(1:12),3] <- pip[2:13,1]*sign(coef.model$postmean[2:13])

6.3.2.3 Prediction plot for marginal effect of conductivity on median infection intensity with Trichodina

newdata1 <- newdata; newdata1[,"Con_av"] <- environment2$Con_av
BMA_Tri_Con_av <- predict(bas.model, newdata = newdata1, estimator = "BMA", se.fit=TRUE)
figure3h = ggplot(environment2, aes(Con_av, BMA_Tri_Con_av$fit)) +
  theme_bw() +
  geom_line(color="red", size=1) +
  geom_ribbon(aes(ymin = (BMA_Tri_Con_av$fit-BMA_Tri_Con_av$se.bma.fit), ymax = (BMA_Tri_Con_av$fit+BMA_Tri_Con_av$se.bma.fit)), alpha = .1) +
  geom_point(data = environment2, aes(x=Con_av, y=medin$Tri)) +
  labs(x=expression(paste("Conductivity [", mu, "S/cm]")), y=expression(paste(italic("Trichodina"), " sp. - infection intensity"))) +
  theme(axis.title.x = element_text(size=12),
        axis.title.y = element_text(size=12)) +  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  labs(subtitle = "h)")
figure3h

6.3.2.3 Prediction plot for marginal effect of COD on median infection intensity with Trichodina

newdata1 <- newdata; newdata1[,"COD_av"] <- environment2$COD_av
BMA_Tri_COD_av <- predict(bas.model, newdata = newdata1, estimator = "BMA", se.fit=TRUE)
figure3d = ggplot(environment2, aes(COD_av, BMA_Tri_COD_av$fit)) +
  theme_bw() +
  geom_line(color="red", size=1) +
  geom_ribbon(aes(ymin = (BMA_Tri_COD_av$fit-BMA_Tri_COD_av$se.bma.fit), ymax = (BMA_Tri_COD_av$fit+BMA_Tri_COD_av$se.bma.fit)), alpha = .1) +
  geom_point(data = environment2, aes(x=COD_av, y=medin$Tri)) +
  labs(x=expression(paste("Chemical Oxygen Demand [mg/L]")), y=expression(paste(italic("Trichodina"), " sp. - infection intensity"))) +
  theme(axis.title.x = element_text(size=12),
        axis.title.y = element_text(size=12)) +  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  labs(subtitle = "d)")
figure3d

6.4.1 Prevalence

bas.model <- bas.lm(prev$Tri ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av 
                     + NH4._av + Nt_av + pool_riffle + meander + netcen + 
                       updist, data=environment2, prior="JZS")
yhat = fitted(bas.model, estimator = "BMA") #these are the fitted values under BMA
r = bas.model$Y - yhat #these are the model residuals
plot(bas.model)

summary(bas.model)
##              P(B != 0 | Y)   model 1   model 2    model 3   model 4   model 5
## Intercept       1.00000000 1.0000000 1.0000000  1.0000000 1.0000000  1.000000
## avlength        0.04005815 0.0000000 0.0000000  0.0000000 0.0000000  0.000000
## avcondition     0.03389439 0.0000000 0.0000000  0.0000000 0.0000000  0.000000
## T_av            0.03654853 0.0000000 0.0000000  0.0000000 0.0000000  0.000000
## O2_sat_av       0.03749123 0.0000000 0.0000000  0.0000000 0.0000000  0.000000
## Con_av          0.20573471 0.0000000 1.0000000  0.0000000 1.0000000  0.000000
## COD_av          0.04629842 0.0000000 0.0000000  0.0000000 0.0000000  0.000000
## NH4._av         0.06690739 0.0000000 0.0000000  1.0000000 0.0000000  0.000000
## Nt_av           0.04604656 0.0000000 0.0000000  0.0000000 0.0000000  1.000000
## pool_riffle1    0.03411173 0.0000000 0.0000000  0.0000000 0.0000000  0.000000
## meander1        0.03657952 0.0000000 0.0000000  0.0000000 0.0000000  0.000000
## netcen          0.07467299 0.0000000 0.0000000  0.0000000 1.0000000  0.000000
## updist          0.07352673 0.0000000 0.0000000  0.0000000 0.0000000  0.000000
## BF                      NA 0.5616705 0.7274201  0.2732671 1.0000000  0.155205
## PostProbs               NA 0.6354000 0.0686000  0.0258000 0.0171000  0.014600
## R2                      NA 0.0000000 0.1264000  0.0750000 0.2249000  0.044000
## dim                     NA 1.0000000 2.0000000  2.0000000 3.0000000  2.000000
## logmarg                 NA 0.0000000 0.2585887 -0.7204656 0.5768398 -1.286169
image(bas.model, rotate=F)

coef.model <- coef(bas.model)
abs(coef.model$postmean)-2*coef.model$postsd > 0
##  [1]  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [13] FALSE
confint(coef.model)
##                       2.5%        97.5%          beta
## Intercept     7.820697e-01 0.8851797650  8.340580e-01
## avlength      0.000000e+00 0.0000000000 -1.453558e-04
## avcondition   0.000000e+00 0.0000000000 -2.610917e-03
## T_av          0.000000e+00 0.0000000000 -1.135202e-05
## O2_sat_av     0.000000e+00 0.0000000000 -1.893466e-05
## Con_av        0.000000e+00 0.0004278091  6.570599e-05
## COD_av        0.000000e+00 0.0000000000  9.451636e-05
## NH4._av       0.000000e+00 0.0191693460  1.713706e-03
## Nt_av        -1.240461e-04 0.0000000000  1.739937e-04
## pool_riffle1  0.000000e+00 0.0000000000  2.104077e-04
## meander1      0.000000e+00 0.0000000000 -6.652944e-04
## netcen       -3.878705e-06 0.0000000000 -3.317353e-07
## updist       -1.482392e-06 0.0000000000 -1.309710e-07
## attr(,"Probability")
## [1] 0.95
## attr(,"class")
## [1] "confint.bas"
plot(confint(coef.model, parm = 2:11))

## NULL
confint <- confint(coef.model, parm = 2:11)
write.table(confint, 'GyroAA.txt', sep="\t")
pip <- summary(bas.model)
PIP[c(1:12),3] <- pip[2:13,1]*sign(coef.model$postmean[2:13])

6.5 Variation in Glugea infection

bas.model <- bas.glm(pa$Glu ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av + 
                       NH4._av + Nt_av + SM_av + pool_riffle + meander + 
                       spavar$netcen + spavar$updist, data=environment2, betaprior=g.prior(100), family=binomial)
yhat = fitted(bas.model, estimator = "BMA") #these are the fitted values under BMA
r = bas.model$Y - yhat #these are the model residuals
plot(bas.model)

summary(bas.model)
##               P(B != 0 | Y)   model 1    model 2   model 3    model 4
## Intercept         1.0000000  1.000000  1.0000000  1.000000  1.0000000
## avlength          0.6661011  0.000000  1.0000000  0.000000  0.0000000
## avcondition       0.4341431  0.000000  1.0000000  1.000000  1.0000000
## T_av              0.2906006  0.000000  0.0000000  0.000000  0.0000000
## O2_sat_av         0.4036499  0.000000  0.0000000  0.000000  0.0000000
## Con_av            0.9885742  1.000000  1.0000000  1.000000  1.0000000
## COD_av            0.5252563  1.000000  1.0000000  0.000000  1.0000000
## NH4._av           0.8087402  1.000000  1.0000000  1.000000  1.0000000
## Nt_av             0.5433594  1.000000  0.0000000  1.000000  1.0000000
## SM_av             0.4432983  1.000000  0.0000000  1.000000  0.0000000
## pool_riffle1      0.4388672  0.000000  0.0000000  0.000000  0.0000000
## meander1          0.9899902  1.000000  1.0000000  1.000000  1.0000000
## spavar$netcen     0.4765747  0.000000  0.0000000  1.000000  0.0000000
## spavar$updist     0.8123047  1.000000  1.0000000  0.000000  1.0000000
## BF                       NA  1.000000  0.8189871  0.769804  0.7351301
## PostProbs                NA  0.046100  0.0397000  0.038800  0.0341000
## R2                       NA  1.000000  1.0000000  1.000000  1.0000000
## dim                      NA  8.000000  8.0000000  8.000000  8.0000000
## logmarg                  NA -5.808403 -6.0080902 -6.070023 -6.1161110
##                  model 5
## Intercept      1.0000000
## avlength       1.0000000
## avcondition    0.0000000
## T_av           0.0000000
## O2_sat_av      1.0000000
## Con_av         1.0000000
## COD_av         0.0000000
## NH4._av        0.0000000
## Nt_av          0.0000000
## SM_av          0.0000000
## pool_riffle1   1.0000000
## meander1       1.0000000
## spavar$netcen  1.0000000
## spavar$updist  1.0000000
## BF             0.5830791
## PostProbs      0.0313000
## R2             1.0000000
## dim            8.0000000
## logmarg       -6.3478357
image(bas.model, rotate=F)

coef.model <- coef(bas.model)
abs(coef.model$postmean)-2*coef.model$postsd > 0
##  [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [13] FALSE FALSE
confint(coef.model)
##                        2.5%        97.5%          beta
## Intercept     -2.902737e+06 2.926170e+06 -6.339692e+03
## avlength      -4.505021e+04 4.199676e+04  9.029383e+01
## avcondition   -9.023933e+05 8.924758e+05  1.496115e+01
## T_av          -2.790922e+04 3.377942e+04  1.699320e+01
## O2_sat_av     -5.691607e+03 6.068530e+03  9.150769e+00
## Con_av        -1.608333e+03 1.623631e+03  4.293194e+00
## COD_av        -1.527908e+04 1.311479e+04  9.080887e+00
## NH4._av       -2.332931e+05 2.150607e+05 -2.537950e+02
## Nt_av         -4.507645e+04 4.741178e+04  3.593276e+01
## SM_av         -2.247258e+03 2.416569e+03 -1.167050e+00
## pool_riffle1  -1.421694e+05 1.397596e+05  8.653441e+01
## meander1      -4.926730e+05 5.194216e+05 -1.344934e+03
## spavar$netcen -1.058765e+01 9.300020e+00 -8.848165e-03
## spavar$updist -4.812767e+00 5.076095e+00 -7.824641e-03
## attr(,"Probability")
## [1] 0.95
## attr(,"class")
## [1] "confint.bas"
plot(confint(coef.model, parm = 2:11))

## NULL
confint <- confint(coef.model, parm = 2:11)
pip <- summary(bas.model)
PIP[c(1:12),9] <- pip[2:13,1]*sign(coef.model$postmean[2:13])

6.5 Variation in Contracaecum infection

bas.model <- bas.glm(pa$Con ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av + 
                       NH4._av + Nt_av + SM_av + pool_riffle + meander + 
                       spavar$netcen + spavar$updist, data=environment2, betaprior=g.prior(100), family=binomial)
yhat = fitted(bas.model, estimator = "BMA") #these are the fitted values under BMA
r = bas.model$Y - yhat #these are the model residuals
plot(bas.model)

summary(bas.model)
##               P(B != 0 | Y)   model 1     model 2     model 3     model 4
## Intercept        1.00000000   1.00000   1.0000000   1.0000000   1.0000000
## avlength         0.08469238   0.00000   1.0000000   0.0000000   0.0000000
## avcondition      0.02604980   0.00000   0.0000000   1.0000000   0.0000000
## T_av             0.01040039   0.00000   0.0000000   0.0000000   0.0000000
## O2_sat_av        0.02716064   0.00000   0.0000000   0.0000000   1.0000000
## Con_av           0.01129150   0.00000   0.0000000   0.0000000   0.0000000
## COD_av           0.01397705   0.00000   0.0000000   0.0000000   0.0000000
## NH4._av          0.01912842   0.00000   0.0000000   0.0000000   0.0000000
## Nt_av            0.01209717   0.00000   0.0000000   0.0000000   0.0000000
## SM_av            0.01187744   0.00000   0.0000000   0.0000000   0.0000000
## pool_riffle1     0.01187744   0.00000   0.0000000   0.0000000   0.0000000
## meander1         0.01021729   0.00000   0.0000000   0.0000000   0.0000000
## spavar$netcen    0.01888428   0.00000   0.0000000   0.0000000   0.0000000
## spavar$updist    0.01214600   0.00000   0.0000000   0.0000000   0.0000000
## BF                       NA   1.00000   0.9507459   0.2581488   0.2583995
## PostProbs                NA   0.78470   0.0608000   0.0152000   0.0145000
## R2                       NA   0.00000   0.0864000   0.0365000   0.0366000
## dim                      NA   1.00000   2.0000000   2.0000000   2.0000000
## logmarg                  NA -25.82594 -25.8764468 -27.1801577 -27.1791867
##                 model 5
## Intercept       1.00000
## avlength        0.00000
## avcondition     0.00000
## T_av            0.00000
## O2_sat_av       0.00000
## Con_av          0.00000
## COD_av          0.00000
## NH4._av         1.00000
## Nt_av           0.00000
## SM_av           0.00000
## pool_riffle1    0.00000
## meander1        0.00000
## spavar$netcen   0.00000
## spavar$updist   0.00000
## BF              0.19791
## PostProbs       0.01250
## R2              0.02640
## dim             2.00000
## logmarg       -27.44588
image(bas.model, rotate=F)

coef.model <- coef(bas.model)
abs(coef.model$postmean)-2*coef.model$postsd > 0
##  [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [13] FALSE FALSE
confint(coef.model)
##                    2.5%     97.5%          beta
## Intercept     -9.938373 2.1628712 -5.683824e-01
## avlength       0.000000 0.1692877  1.573802e-02
## avcondition    0.000000 0.0000000 -1.691605e-01
## T_av           0.000000 0.0000000  7.921083e-04
## O2_sat_av      0.000000 0.0000000  1.190308e-03
## Con_av         0.000000 0.0000000  5.101715e-06
## COD_av         0.000000 0.0000000 -3.712753e-04
## NH4._av        0.000000 0.0000000 -5.964150e-03
## Nt_av          0.000000 0.0000000 -9.152172e-04
## SM_av          0.000000 0.0000000  1.162131e-04
## pool_riffle1   0.000000 0.0000000 -5.259321e-03
## meander1       0.000000 0.0000000  2.907684e-03
## spavar$netcen  0.000000 0.0000000 -8.702274e-07
## spavar$updist  0.000000 0.0000000 -1.184043e-07
## attr(,"Probability")
## [1] 0.95
## attr(,"class")
## [1] "confint.bas"
plot(confint(coef.model, parm = 2:11))

## NULL
confint <- confint(coef.model, parm = 2:11)
pip <- summary(bas.model)
PIP[c(1:12),10] <- pip[2:13,1]*sign(coef.model$postmean[2:13])

Variation in Anguillicoloides infection

bas.model <- bas.glm(pa$Ang ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av + 
                       NH4._av + Nt_av + SM_av + pool_riffle + meander + 
                       spavar$netcen + spavar$updist, data=environment2, betaprior=g.prior(100), family=binomial)
yhat = fitted(bas.model, estimator = "BMA") #these are the fitted values under BMA
r = bas.model$Y - yhat #these are the model residuals
plot(bas.model)

summary(bas.model)
##               P(B != 0 | Y)   model 1     model 2     model 3     model 4
## Intercept         1.0000000  1.000000   1.0000000   1.0000000   1.0000000
## avlength          0.9313965  1.000000   1.0000000   1.0000000   1.0000000
## avcondition       0.3154663  0.000000   0.0000000   0.0000000   0.0000000
## T_av              0.9406494  1.000000   1.0000000   1.0000000   1.0000000
## O2_sat_av         0.8143555  1.000000   1.0000000   1.0000000   1.0000000
## Con_av            0.5466187  0.000000   1.0000000   1.0000000   0.0000000
## COD_av            0.3676270  0.000000   0.0000000   1.0000000   0.0000000
## NH4._av           0.9558350  1.000000   1.0000000   1.0000000   1.0000000
## Nt_av             0.5304565  0.000000   1.0000000   0.0000000   0.0000000
## SM_av             0.9200195  1.000000   1.0000000   1.0000000   1.0000000
## pool_riffle1      0.3712646  0.000000   0.0000000   0.0000000   1.0000000
## meander1          0.9378540  1.000000   1.0000000   1.0000000   1.0000000
## spavar$netcen     0.8593750  1.000000   1.0000000   1.0000000   1.0000000
## spavar$updist     0.5289795  1.000000   0.0000000   0.0000000   1.0000000
## BF                       NA  1.000000   0.4728289   0.3794633   0.2858158
## PostProbs                NA  0.079300   0.0720000   0.0531000   0.0442000
## R2                       NA  1.000000   1.0000000   1.0000000   1.0000000
## dim                      NA  9.000000  10.0000000  10.0000000  10.0000000
## logmarg                  NA -9.802499 -10.5515211 -10.7714968 -11.0549070
##                   model 5
## Intercept       1.0000000
## avlength        1.0000000
## avcondition     1.0000000
## T_av            1.0000000
## O2_sat_av       0.0000000
## Con_av          1.0000000
## COD_av          0.0000000
## NH4._av         1.0000000
## Nt_av           1.0000000
## SM_av           1.0000000
## pool_riffle1    0.0000000
## meander1        1.0000000
## spavar$netcen   1.0000000
## spavar$updist   0.0000000
## BF              0.2037316
## PostProbs       0.0299000
## R2              1.0000000
## dim            10.0000000
## logmarg       -11.3934513
image(bas.model, rotate=F)

coef.model <- coef(bas.model)
abs(coef.model$postmean)-2*coef.model$postsd > 0
##  [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [13] FALSE FALSE
confint(coef.model)
##                        2.5%        97.5%          beta
## Intercept     -2.493583e+06 2.220172e+06 -6.314576e+03
## avlength      -4.633181e+04 4.189389e+04  1.289938e+02
## avcondition   -8.259937e+05 7.290298e+05 -2.425364e+02
## T_av          -1.617354e+05 1.522072e+05  5.783392e+02
## O2_sat_av     -1.894797e+04 1.783841e+04 -3.694674e+01
## Con_av        -6.725601e+02 6.298949e+02  3.253851e-01
## COD_av        -7.695858e+03 7.302203e+03  1.422581e+00
## NH4._av       -3.522546e+05 3.350507e+05 -1.031216e+03
## Nt_av         -4.309152e+04 4.253718e+04 -4.515461e+01
## SM_av         -5.880688e+03 5.728770e+03  1.263286e+01
## pool_riffle1  -1.148759e+05 1.313316e+05 -1.031862e+02
## meander1      -3.968940e+05 3.849816e+05 -1.367302e+03
## spavar$netcen -1.525973e+01 1.372595e+01 -2.731076e-02
## spavar$updist -3.007013e+00 2.830525e+00  1.319948e-03
## attr(,"Probability")
## [1] 0.95
## attr(,"class")
## [1] "confint.bas"
plot(confint(coef.model, parm = 2:11))

## NULL
confint <- confint(coef.model, parm = 2:11)
pip <- summary(bas.model)
PIP[c(1:12),11] <- pip[2:13,1]*sign(coef.model$postmean[2:13])

Variation in Individual Parasitation Index (all parasites)

bas.model <- bas.lm(avPI ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av 
                     + NH4._av + Nt_av + pool_riffle + meander + netcen + 
                       updist, data=environment2, prior="JZS")
yhat = fitted(bas.model, estimator = "BMA") #these are the fitted values under BMA
r = bas.model$Y - yhat #these are the model residuals
plot(bas.model)

summary(bas.model)
##              P(B != 0 | Y)    model 1    model 2    model 3  model 4  model 5
## Intercept        1.0000000 1.00000000 1.00000000 1.00000000 1.000000 1.000000
## avlength         0.2261960 0.00000000 0.00000000 0.00000000 0.000000 0.000000
## avcondition      0.1810890 0.00000000 0.00000000 0.00000000 0.000000 0.000000
## T_av             0.4574073 0.00000000 1.00000000 0.00000000 1.000000 0.000000
## O2_sat_av        0.1952767 0.00000000 0.00000000 0.00000000 0.000000 0.000000
## Con_av           0.4598919 0.00000000 0.00000000 0.00000000 0.000000 1.000000
## COD_av           0.5684747 0.00000000 0.00000000 0.00000000 0.000000 1.000000
## NH4._av          0.2748144 0.00000000 0.00000000 0.00000000 0.000000 0.000000
## Nt_av            0.4582802 0.00000000 0.00000000 1.00000000 1.000000 0.000000
## pool_riffle1     0.3054914 0.00000000 0.00000000 0.00000000 0.000000 0.000000
## meander1         0.1985617 0.00000000 0.00000000 0.00000000 0.000000 0.000000
## netcen           0.1910872 0.00000000 0.00000000 0.00000000 0.000000 0.000000
## updist           0.4153657 0.00000000 0.00000000 0.00000000 0.000000 1.000000
## BF                      NA 0.01969453 0.09251289 0.08737546 0.355023 1.000000
## PostProbs               NA 0.10850000 0.04250000 0.04010000 0.029600 0.025000
## R2                      NA 0.00000000 0.18980000 0.18710000 0.325400 0.424800
## dim                     NA 1.00000000 2.00000000 2.00000000 3.000000 4.000000
## logmarg                 NA 0.00000000 1.54700725 1.48987379 2.891842 3.927415
image(bas.model, rotate=F)

coef.model <- coef(bas.model)
abs(coef.model$postmean)-2*coef.model$postsd > 0
##  [1]  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [13] FALSE
confint(coef.model)
##                       2.5%        97.5%          beta
## Intercept     1.139112e+00 1.732523e+00  1.422090e+00
## avlength     -8.655588e-02 1.877262e-02 -7.992245e-03
## avcondition  -2.976576e+00 1.529684e+00 -9.325811e-02
## T_av         -7.038395e-04 3.716692e-01  8.644408e-02
## O2_sat_av    -2.345726e-02 1.521722e-02 -7.476618e-04
## Con_av       -3.164277e-06 3.173245e-03  8.405554e-04
## COD_av        0.000000e+00 6.148546e-02  2.109188e-02
## NH4._av      -5.306465e-01 1.492837e-01 -4.965441e-02
## Nt_av         0.000000e+00 2.234194e-01  5.507750e-02
## pool_riffle1 -3.737605e-02 8.953177e-01  1.389140e-01
## meander1     -5.590964e-01 2.250265e-01 -3.866036e-02
## netcen       -2.998202e-05 1.708883e-05 -1.270941e-06
## updist       -2.627701e-05 0.000000e+00 -5.952452e-06
## attr(,"Probability")
## [1] 0.95
## attr(,"class")
## [1] "confint.bas"
plot(confint(coef.model, parm = 2:11))

## NULL
confint <- confint(coef.model, parm = 2:11)
pip <- summary(bas.model)
PIP[c(1:12),12] <- pip[2:13,1]*sign(coef.model$postmean[2:13])

Prediction plot for marginal effect of COD on Individual Parasitation Index

newdata1 <- newdata; newdata1[,"COD_av"] <- environment2$COD_av
BMA_IPI_COD_av <- predict(bas.model, newdata = newdata1, estimator = "BMA", se.fit=TRUE)
figure3b = ggplot(environment2, aes(COD_av, BMA_IPI_COD_av$fit)) +
  theme_bw() +
  geom_line(color="red", size=1) +
  geom_ribbon(aes(ymin = (BMA_IPI_COD_av$fit-BMA_IPI_COD_av$se.bma.fit), ymax = (BMA_IPI_COD_av$fit+BMA_IPI_COD_av$se.bma.fit)), alpha = .1) +
  geom_point(data = environment2, aes(x=COD_av, y=avPI)) +
  labs(x=expression("Chemical Oxygen Demand [mg/L]"), y=expression("I"[PI]*" - all parasites")) +
  theme(axis.title.x = element_text(size=12),
        axis.title.y = element_text(size=12)) +  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  labs(subtitle = "b)")
figure3b

Variation in Individual Parasitation Index (only ectoparasites)

bas.model <- bas.lm(avPI_ecto ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av 
                     + NH4._av + Nt_av + pool_riffle + meander + netcen + 
                       updist, data=environment2, prior="JZS")
yhat = fitted(bas.model, estimator = "BMA") #these are the fitted values under BMA
r = bas.model$Y - yhat #these are the model residuals
plot(bas.model)

summary(bas.model)
##              P(B != 0 | Y)   model 1      model 2   model 3  model 4  model 5
## Intercept        1.0000000 1.0000000  1.000000000 1.0000000 1.000000 1.000000
## avlength         0.4669536 0.0000000  1.000000000 0.0000000 0.000000 1.000000
## avcondition      0.3758074 0.0000000  1.000000000 0.0000000 0.000000 0.000000
## T_av             0.2753561 0.0000000  1.000000000 0.0000000 0.000000 0.000000
## O2_sat_av        0.2919614 0.0000000  1.000000000 0.0000000 0.000000 0.000000
## Con_av           0.6617062 0.0000000  1.000000000 1.0000000 0.000000 0.000000
## COD_av           0.8437205 0.0000000  1.000000000 1.0000000 1.000000 1.000000
## NH4._av          0.4163705 0.0000000  1.000000000 0.0000000 0.000000 0.000000
## Nt_av            0.8302866 1.0000000  1.000000000 0.0000000 1.000000 1.000000
## pool_riffle1     0.4025108 0.0000000  1.000000000 0.0000000 0.000000 0.000000
## meander1         0.5031625 0.0000000  1.000000000 0.0000000 0.000000 0.000000
## netcen           0.3010717 0.0000000  1.000000000 0.0000000 0.000000 0.000000
## updist           0.3285744 0.0000000  1.000000000 0.0000000 0.000000 0.000000
## BF                      NA 0.1318791  0.008313608 0.4282678 0.319205 1.000000
## PostProbs               NA 0.0510000  0.038600000 0.0301000 0.022500 0.021100
## R2                      NA 0.3007000  0.682800000 0.4141000 0.403700 0.496600
## dim                     NA 2.0000000 13.000000000 3.0000000 3.000000 4.000000
## logmarg                 NA 4.0635523  1.299560343 5.2414154 4.947500 6.089422
image(bas.model, rotate=F)

coef.model <- coef(bas.model)
abs(coef.model$postmean)-2*coef.model$postsd > 0
##  [1]  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [13] FALSE
confint(coef.model)
##                       2.5%        97.5%          beta
## Intercept     6.885848e-01 1.160730e+00  9.199577e-01
## avlength     -1.053581e-01 1.014527e-02 -2.357999e-02
## avcondition  -4.323979e+00 6.216264e-01 -7.055549e-01
## T_av         -1.006276e-01 1.379698e-01  5.818763e-03
## O2_sat_av    -1.303563e-02 2.535981e-02  1.122707e-03
## Con_av       -2.688244e-06 3.154139e-03  1.203349e-03
## COD_av        0.000000e+00 6.418310e-02  3.405818e-02
## NH4._av      -4.708824e-01 1.076874e-01 -8.034086e-02
## Nt_av         0.000000e+00 2.533901e-01  1.257830e-01
## pool_riffle1 -8.969849e-02 7.899060e-01  1.490513e-01
## meander1     -9.360961e-01 1.708101e-02 -2.358742e-01
## netcen       -3.378862e-05 2.125678e-05 -2.447510e-06
## updist       -1.563006e-05 4.018463e-06 -1.836166e-06
## attr(,"Probability")
## [1] 0.95
## attr(,"class")
## [1] "confint.bas"
plot(confint(coef.model, parm = 2:11))

## NULL
confint <- confint(coef.model, parm = 2:11)
pip <- summary(bas.model)
PIP[c(1:12),13] <- pip[2:13,1]*sign(coef.model$postmean[2:13])

Prediction plot for marginal effect of COD on Individual Parasitation Index (only ectoparasites)

newdata1 <- newdata; newdata1[,"COD_av"] <- environment2$COD_av
BMA_IPIecto_COD_av <- predict(bas.model, newdata = newdata1, estimator = "BMA", se.fit=TRUE)
figure3c = ggplot(environment2, aes(COD_av, BMA_IPIecto_COD_av$fit)) +
  theme_bw() +
  geom_line(color="red", size=1) +
  geom_ribbon(aes(ymin = (BMA_IPIecto_COD_av$fit-BMA_IPIecto_COD_av$se.bma.fit), ymax = (BMA_IPIecto_COD_av$fit+BMA_IPIecto_COD_av$se.bma.fit)), alpha = .1) +
  geom_point(data = environment2, aes(x=COD_av, y=avPI_ecto)) +
  labs(x=expression("Chemical Oxygen Demand [mg/L]"), y=expression("I"[PI]*" - ectoparasites")) +
  theme(axis.title.x = element_text(size=12),
        axis.title.y = element_text(size=12)) +  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  labs(subtitle = "c)")
figure3c

Prediction plot for marginal effect of total nitrogen on Individual Parasitation Index (only ectoparasites)

newdata1 <- newdata; newdata1[,"Nt_av"] <- environment2$Nt_av
BMA_IPIecto_Nt_av <- predict(bas.model, newdata = newdata1, estimator = "BMA", se.fit=TRUE)
figure3f = ggplot(environment2, aes(Nt_av, BMA_IPIecto_Nt_av$fit)) +
  theme_bw() +
  geom_line(color="red", size=1) +
  geom_ribbon(aes(ymin = (BMA_IPIecto_Nt_av$fit-BMA_IPIecto_Nt_av$se.bma.fit), ymax = (BMA_IPIecto_Nt_av$fit+BMA_IPIecto_Nt_av$se.bma.fit)), alpha = .1) +
  geom_point(data = environment2, aes(x=Nt_av, y=avPI_ecto)) +
  labs(x=expression("Nt [mg/L]"), y=expression("I"[PI]*" - ectoparasites")) +
  theme(axis.title.x = element_text(size=12),
        axis.title.y = element_text(size=12)) +  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  labs(subtitle = "f)")
figure3f

Prediction plot for marginal effect of conductivity on Individual Parasitation Index (only ectoparasites)

newdata1 <- newdata; newdata1[,"Con_av"] <- environment2$Con_av
BMA_IPIecto_Con_av <- predict(bas.model, newdata = newdata1, estimator = "BMA", se.fit=TRUE)
figure3g = ggplot(environment2, aes(Con_av, BMA_IPIecto_Con_av$fit)) +
  theme_bw() +
  geom_line(color="red", size=1) +
  geom_ribbon(aes(ymin = (BMA_IPIecto_Con_av$fit-BMA_IPIecto_Con_av$se.bma.fit), ymax = (BMA_IPIecto_Con_av$fit+BMA_IPIecto_Con_av$se.bma.fit)), alpha = .1) +
  geom_point(data = environment2, aes(x=Con_av, y=avPI)) +
  labs(x=expression(paste("Conductivity [", mu, "S/cm]")), y=expression("I"[PI]*" - all parasites")) +
  theme(axis.title.x = element_text(size=12),
        axis.title.y = element_text(size=12)) +  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  labs(subtitle = "g)")
figure3g

Prediction plot for marginal effect of meanders on Individual Parasitation Index (only ectoparasites)

newdata1 <- newdata; newdata1[,"meander"] <- environment2$meander
BMA_IPIecto_meander <- predict(bas.model, newdata = newdata1, estimator = "BMA", se.fit=TRUE)
figure3l = ggplot(environment2, aes(meander, BMA_IPIecto_meander$fit)) +
  theme_bw() +
  #geom_line(color="red", size=1) +
  #geom_ribbon(aes(ymin = (BMA$fit-BMA$se.bma.fit), ymax = (BMA$fit+BMA$se.bma.fit)), alpha = .1) +
  geom_point(data = environment2, aes(x=meander, y=avPI)) +
  geom_boxplot(aes(lower = (BMA_IPIecto_meander$fit-BMA_IPIecto_meander$se.bma.fit), middle = BMA_IPIecto_meander$fit, upper = (BMA_IPIecto_meander$fit+BMA_IPIecto_meander$se.bma.fit))) +
  labs(x=expression("Meander"), y=expression("I"[PI]*" - ectoparasites")) +
  theme(axis.title.x = element_text(size=12),
        axis.title.y = element_text(size=12)) +  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  labs(subtitle = "l)")
figure3l

Variation in Individual Parasitation Index (only endoparasites)

bas.model <- bas.lm(avPI_endo ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av 
                     + NH4._av + Nt_av + pool_riffle + meander + netcen + 
                       updist, data=environment2, prior="JZS")
yhat = fitted(bas.model, estimator = "BMA") #these are the fitted values under BMA
r = bas.model$Y - yhat #these are the model residuals
plot(bas.model)

summary(bas.model)
##              P(B != 0 | Y)   model 1  model 2     model 3    model 4   model 5
## Intercept       1.00000000 1.0000000 1.000000  1.00000000  1.0000000 1.0000000
## avlength        0.06442140 0.0000000 0.000000  0.00000000  0.0000000 0.0000000
## avcondition     0.12246410 0.0000000 0.000000  1.00000000  0.0000000 0.0000000
## T_av            0.30397004 0.0000000 1.000000  0.00000000  0.0000000 1.0000000
## O2_sat_av       0.05178919 0.0000000 0.000000  0.00000000  0.0000000 0.0000000
## Con_av          0.05158494 0.0000000 0.000000  0.00000000  0.0000000 0.0000000
## COD_av          0.04593088 0.0000000 0.000000  0.00000000  0.0000000 0.0000000
## NH4._av         0.08528364 0.0000000 0.000000  0.00000000  0.0000000 1.0000000
## Nt_av           0.06477759 0.0000000 0.000000  0.00000000  0.0000000 0.0000000
## pool_riffle1    0.04504579 0.0000000 0.000000  0.00000000  0.0000000 0.0000000
## meander1        0.05899866 0.0000000 0.000000  0.00000000  0.0000000 0.0000000
## netcen          0.05414430 0.0000000 0.000000  0.00000000  0.0000000 0.0000000
## updist          0.10936917 0.0000000 0.000000  0.00000000  1.0000000 0.0000000
## BF                      NA 0.3484238 1.000000  0.33059983  0.2877395 0.6134470
## PostProbs               NA 0.4921000 0.117700  0.03890000  0.0339000 0.0131000
## R2                      NA 0.0000000 0.166100  0.11040000  0.1032000 0.2244000
## dim                     NA 1.0000000 2.000000  2.00000000  2.0000000 3.0000000
## logmarg                 NA 0.0000000 1.054336 -0.05251095 -0.1913639 0.5656742
image(bas.model, rotate=F)

coef.model <- coef(bas.model)
abs(coef.model$postmean)-2*coef.model$postsd > 0
##  [1]  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [13] FALSE
confint(coef.model)
##                       2.5%        97.5%          beta
## Intercept     3.062561e-01 6.936441e-01  4.947994e-01
## avlength     -1.350725e-04 1.781226e-02  1.398963e-03
## avcondition   0.000000e+00 2.576839e+00  2.446911e-01
## T_av          0.000000e+00 2.160719e-01  4.520999e-02
## O2_sat_av    -3.667475e-04 1.264879e-04 -2.446396e-04
## Con_av        0.000000e+00 0.000000e+00 -1.385311e-05
## COD_av        0.000000e+00 0.000000e+00 -1.072193e-04
## NH4._av      -1.139164e-01 0.000000e+00 -9.322620e-03
## Nt_av        -2.667990e-02 2.022951e-04 -2.149910e-03
## pool_riffle1  0.000000e+00 0.000000e+00  2.267863e-03
## meander1      0.000000e+00 8.606124e-02  1.005422e-02
## netcen       -1.545415e-07 1.173066e-06 -3.607730e-07
## updist       -8.998181e-06 0.000000e+00 -8.091451e-07
## attr(,"Probability")
## [1] 0.95
## attr(,"class")
## [1] "confint.bas"
plot(confint(coef.model, parm = 2:11))
## Warning in arrows(x[not.deg], ci[not.deg, 1], x[not.deg], ci[not.deg, 2], :
## zero-length arrow is of indeterminate angle and so skipped

## NULL
confint <- confint(coef.model, parm = 2:11)
pip <- summary(bas.model)
PIP[c(1:12),14] <- pip[2:13,1]*sign(coef.model$postmean[2:13])

Figure 3

library(gridExtra)
pdf("Figure3.pdf", width = 10.8, height = 14.4)
grid.arrange(figure3a, figure3b, figure3c, figure3d, figure3e, figure3f, figure3g, figure3h, figure3i, figure3j, figure3k, figure3l)
dev.off()
## png 
##   2

7. BORAL analysis

Model-based analysis of multivariate abundance data using Bayesian Markov chain Monte Carlo methods for parameter estimation

7.1 BORAL analysis for average abundances of parasites

data$Site <- as.factor(data$site)
levels(data$site) <- levels(as.factor(environment2$Site))
data_m <- merge(data, environment2, by = "Site")
data_all <- na.omit(data_m) 
names(data_all)
##  [1] "Site"                "site"                "fish"               
##  [4] "parspeciesrichness"  "div_shannon"         "div_simpson"        
##  [7] "temp"                "pH"                  "conductivity"       
## [10] "nitrogen"            "phosphorus"          "oxygen"             
## [13] "netcen.x"            "updist.x"            "updist2"            
## [16] "updist3"             "fishspeciesrichness" "weight"             
## [19] "weigh..g."           "length"              "SMI"                
## [22] "Sex"                 "Gyr"                 "Tri"                
## [25] "Glu"                 "ecto_screener"       "Con"                
## [28] "CysL"                "Pro"                 "Aca"                
## [31] "Cam"                 "Ang"                 "CysI"               
## [34] "endo_screener"       "PI"                  "PI_ecto"            
## [37] "PI_endo"             "T_av"                "O2_sat_av"          
## [40] "Con_av"              "COD_av"              "NH4._av"            
## [43] "Nt_av"               "SM_av"               "pool_riffle"        
## [46] "meander"             "updist.y"            "netcen.y"
avcondition <- aggregate(data$SMI, by = list(data[,1]), function(x){mean(x, na.rm =T)})[,2]
avlength <- aggregate(data$length, by = list(data[,1]), function(x){mean(x, na.rm =T)})[,2]

y <- round(cbind(avab$Gyr, avab$Tri, avab$Glu, avab$Con, avab$Ang))
X <- cbind(avcondition,
           avlength,
           environment2$T_av,
           environment2$O2_sat_av,
           environment2$Con_av,
           environment2$COD_av, 
           environment2$NH4._av, 
           environment2$Nt_av, 
           environment2$netcen, 
           environment2$updist, 
           as.numeric(environment2$pool_riffle), 
           as.numeric(environment2$meander))
colnames(X) <- c("avcondition", "avlength", "T", "O2", "Con", "COD", "NH4", "Nt", "netcen", "updist", "pool_riffle", "meander")

example_mcmc_control <- list(n.burnin = 1000, n.iteration = 10000, n.thin = 1)
testpath <- file.path(tempdir(), "jagsboralmodel.txt")
paramod <- boral(y, X = X,
                      family = "negative.binomial",
                      mcmc.control = example_mcmc_control,
                      model.name = testpath,
                      lv.control = list(num.lv = 2, type = "independent"),
                      save.model = TRUE)
## module glm loaded
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 185
##    Unobserved stochastic nodes: 338
##    Total graph size: 2173
## 
## Initializing model
plot(paramod)
## NULL

coefsplot(covname = "avcondition", object = paramod) #Condition

coefsplot(covname = "avlength", object = paramod) #Length

coefsplot(covname = "T", object = paramod) #Temperature

coefsplot(covname = "O2", object = paramod) #Oxygen

coefsplot(covname = "Con", object = paramod) #Conductivity

coefsplot(covname = "COD", object = paramod) #COD

coefsplot(covname = "NH4", object = paramod) #NH4

coefsplot(covname = "Nt", object = paramod) #Nt

coefsplot(covname = "netcen", object = paramod) #netcen

coefsplot(covname = "updist", object = paramod) #updist

coefsplot(covname = "pool_riffle", object = paramod) #poolriffle

coefsplot(covname = "meander", object = paramod) #meander

envcors <- get.enviro.cor(paramod)
rescors <- get.residual.cor(paramod)
corrplot(envcors$sig.cor, type = "lower", diag = FALSE, title = "Correlations due to covariates", mar = c(3,0.5,2,1), tl.srt = 45)

corrplot(rescors$sig.cor, type = "lower", diag = FALSE, title = "Residual correlations", mar = c(3,0.5,2,1), tl.srt = 45)

7.1 BORAL analysis for median infection intensities of parasites

y <- round(cbind(medin$Gyr, medin$Tri, medin$Glu, medin$Con, medin$Ang))
paramod <- boral(y, X = X,
                      family = "negative.binomial",
                      mcmc.control = example_mcmc_control,
                      model.name = testpath,
                      lv.control = list(num.lv = 2, type = "independent"),
                      save.model = TRUE)
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 185
##    Unobserved stochastic nodes: 338
##    Total graph size: 2173
## 
## Initializing model
plot(paramod)
## NULL

coefsplot(covname = "avcondition", object = paramod) #Condition

coefsplot(covname = "avlength", object = paramod) #Length

coefsplot(covname = "T", object = paramod) #Temperature

coefsplot(covname = "O2", object = paramod) #Oxygen

coefsplot(covname = "Con", object = paramod) #Conductivity

coefsplot(covname = "COD", object = paramod) #COD

coefsplot(covname = "NH4", object = paramod) #NH4

coefsplot(covname = "Nt", object = paramod) #Nt

coefsplot(covname = "netcen", object = paramod) #netcen

coefsplot(covname = "updist", object = paramod) #updist

coefsplot(covname = "pool_riffle", object = paramod) #poolriffle

coefsplot(covname = "meander", object = paramod) #meander

envcors <- get.enviro.cor(paramod)
rescors <- get.residual.cor(paramod)
corrplot(envcors$sig.cor, type = "lower", diag = FALSE, title = "Correlations due to covariates", mar = c(3,0.5,2,1), tl.srt = 45)

corrplot(rescors$sig.cor, type = "lower", diag = FALSE, title = "Residual correlations", mar = c(3,0.5,2,1), tl.srt = 45)

8. Multivariate analysis

8.1 Component communities

# Component communities: Bray-Curtis dissimilarities based on Hellinger transformed average abundance data
spe.hel_bray <- vegdist(decostand(avab[,-1], na.rm=T, method="hellinger"), method="bray", na.rm=T)

# Check whether Euclidean and Bray-Curtis distances are comparable
spe.hel_euc <- vegdist(decostand(avab[,-1], na.rm=T, method="hellinger"), method="euc", na.rm=T)
plot(spe.hel_bray, spe.hel_euc)

mantel(spe.hel_bray, spe.hel_euc)
## 
## Mantel statistic based on Pearson's product-moment correlation 
## 
## Call:
## mantel(xdis = spe.hel_bray, ydis = spe.hel_euc) 
## 
## Mantel statistic r: 0.9648 
##       Significance: 0.001 
## 
## Upper quantiles of permutations (null model):
##   90%   95% 97.5%   99% 
## 0.110 0.144 0.186 0.206 
## Permutation: free
## Number of permutations: 999
model_adonis = adonis(spe.hel_bray ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av 
                     + NH4._av + Nt_av + pool_riffle + meander + netcen + 
                       updist, data=environment2)
## 'adonis' will be deprecated: use 'adonis2' instead
model_adonis$aov.tab
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##             Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)  
## avlength     1    0.1863 0.18629  1.8398 0.04303  0.110  
## avcondition  1    0.1305 0.13055  1.2893 0.03016  0.272  
## T_av         1    0.3204 0.32039  3.1643 0.07401  0.020 *
## O2_sat_av    1    0.1368 0.13678  1.3509 0.03160  0.266  
## Con_av       1    0.1326 0.13262  1.3098 0.03064  0.280  
## COD_av       1    0.0657 0.06568  0.6486 0.01517  0.658  
## NH4._av      1    0.2040 0.20399  2.0146 0.04712  0.110  
## Nt_av        1    0.1451 0.14509  1.4329 0.03352  0.212  
## pool_riffle  1    0.0853 0.08529  0.8424 0.01970  0.496  
## meander      1    0.2029 0.20286  2.0035 0.04686  0.098 .
## netcen       1    0.2173 0.21728  2.1459 0.05019  0.071 .
## updist       1    0.0719 0.07193  0.7104 0.01662  0.577  
## Residuals   24    2.4301 0.10125         0.56137         
## Total       36    4.3288                 1.00000         
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# environmental variables
env_select <- environment2[,c("T_av", "O2_sat_av", "Con_av", "COD_av", "NH4._av", "Nt_av", "pool_riffle", "meander", "netcen", "updist")]
env_select$pool_riffle <- as.numeric(env_select$pool_riffle)
env_select$meander <- as.numeric(env_select$meander)

pca <- prcomp(env_select, scale.=T)
summary(pca)
## Importance of components:
##                           PC1    PC2    PC3    PC4     PC5     PC6     PC7
## Standard deviation     1.7124 1.5545 1.1221 1.0140 0.88807 0.79463 0.56647
## Proportion of Variance 0.2933 0.2416 0.1259 0.1028 0.07887 0.06314 0.03209
## Cumulative Proportion  0.2933 0.5349 0.6608 0.7636 0.84248 0.90563 0.93771
##                            PC8     PC9    PC10
## Standard deviation     0.50483 0.46939 0.38429
## Proportion of Variance 0.02549 0.02203 0.01477
## Cumulative Proportion  0.96320 0.98523 1.00000
plot(pca)

biplot(pca)

# Assess the effect of environmental variables on parasite component community dissimilarities using distance based RDA
spe.rda <- dbrda(spe.hel_bray ~ T_av + O2_sat_av + Con_av + COD_av 
                     + NH4._av + Nt_av + pool_riffle + meander, data = env_select)
anova(spe.rda)
## Permutation test for dbrda under reduced model
## Permutation: free
## Number of permutations: 999
## 
## Model: dbrda(formula = spe.hel_bray ~ T_av + O2_sat_av + Con_av + COD_av + NH4._av + Nt_av + pool_riffle + meander, data = env_select)
##          Df SumOfSqs      F Pr(>F)
## Model     8   1.1664 1.2909  0.163
## Residual 28   3.1624
RsquareAdj(spe.rda)$adj.r.squared
## [1] 0.06072755
mod0 <- dbrda(spe.hel_bray ~ 1, env_select[,-c(9:10)])  # Model with intercept only  #edit_PH
mod1 <- dbrda(spe.hel_bray ~ ., env_select[,-c(9:10)])  # Model with all explanatory variables  #edit_PH
step.res <- ordiR2step(mod0, mod1, direction = "both",perm.max = 200)
## Step: R2.adj= 0 
## Call: spe.hel_bray ~ 1 
##  
##                  R2.adjusted
## <All variables>  0.060727547
## + T_av           0.036398792
## + NH4._av        0.020208612
## + meander        0.018502880
## + O2_sat_av      0.004277611
## + Con_av         0.001872668
## + pool_riffle    0.001742860
## <none>           0.000000000
## + Nt_av         -0.002060170
## + COD_av        -0.017968936
## 
##        Df    AIC      F Pr(>F)  
## + T_av  1 54.788 2.3599  0.048 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: R2.adj= 0.03639879 
## Call: spe.hel_bray ~ T_av 
##  
##                 R2.adjusted
## + NH4._av        0.07975512
## <All variables>  0.06072755
## + meander        0.05721602
## + Con_av         0.04792452
## + Nt_av          0.04664053
## + O2_sat_av      0.04174670
## + pool_riffle    0.03938572
## <none>           0.03639879
## + COD_av         0.02324854
step.res$anova  # Summary table
##                   R2.adj Df    AIC      F Pr(>F)  
## + T_av          0.036399  1 54.788 2.3599  0.048 *
## <All variables> 0.060728                          
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(spe.rda, scaling = 1) # it is for technical reasons not possible to plot both site and species scores

summary(spe.rda)
## 
## Call:
## dbrda(formula = spe.hel_bray ~ T_av + O2_sat_av + Con_av + COD_av +      NH4._av + Nt_av + pool_riffle + meander, data = env_select) 
## 
## Partitioning of squared Bray distance:
##               Inertia Proportion
## Total           4.329     1.0000
## Constrained     1.166     0.2695
## Unconstrained   3.162     0.7305
## 
## Eigenvalues, and their contribution to the squared Bray distance 
## 
## Importance of components:
##                       dbRDA1 dbRDA2  dbRDA3  dbRDA4   dbRDA5    dbRDA6
## Eigenvalue            0.5740 0.3372 0.17566 0.09415 0.042525 0.0020670
## Proportion Explained  0.1326 0.0779 0.04058 0.02175 0.009824 0.0004775
## Cumulative Proportion     NA     NA      NA      NA       NA        NA
##                         idbRDA1  idbRDA2   MDS1   MDS2   MDS3    MDS4    MDS5
## Eigenvalue            -0.014797 -0.04441 1.3249 0.8027 0.4680 0.31323 0.29121
## Proportion Explained   0.003418  0.01026 0.3061 0.1854 0.1081 0.07236 0.06727
## Cumulative Proportion        NA       NA     NA     NA     NA      NA      NA
##                          MDS6    MDS7    MDS8    MDS9   MDS10   MDS11   MDS12
## Eigenvalue            0.13644 0.11438 0.08987 0.07065 0.06425 0.02766 0.01576
## Proportion Explained  0.03152 0.02642 0.02076 0.01632 0.01484 0.00639 0.00364
## Cumulative Proportion      NA      NA      NA      NA      NA      NA      NA
##                          MDS13     MDS14      iMDS1     iMDS2    iMDS3
## Eigenvalue            0.011872 2.725e-04 -0.0020314 -0.007356 -0.01078
## Proportion Explained  0.002743 6.295e-05  0.0004693  0.001699  0.00249
## Cumulative Proportion       NA        NA         NA        NA       NA
##                           iMDS4    iMDS5     iMDS6     iMDS7    iMDS8     iMDS9
## Eigenvalue            -0.014744 -0.01935 -0.022495 -0.029849 -0.03342 -0.041489
## Proportion Explained   0.003406  0.00447  0.005197  0.006895  0.00772  0.009584
## Cumulative Proportion        NA       NA        NA        NA       NA        NA
##                         iMDS10   iMDS11   iMDS12   iMDS13   iMDS14
## Eigenvalue            -0.05984 -0.06439 -0.06995 -0.09116 -0.10203
## Proportion Explained   0.01382  0.01487  0.01616  0.02106  0.02357
## Cumulative Proportion       NA       NA       NA       NA       NA
## 
## Accumulated constrained eigenvalues
## Importance of components:
##                       dbRDA1 dbRDA2 dbRDA3  dbRDA4  dbRDA5   dbRDA6  idbRDA1
## Eigenvalue            0.5740 0.3372 0.1757 0.09415 0.04252 0.002067 -0.01480
## Proportion Explained  0.4921 0.2891 0.1506 0.08072 0.03646 0.001772  0.01269
## Cumulative Proportion     NA     NA     NA      NA      NA       NA       NA
##                        idbRDA2
## Eigenvalue            -0.04441
## Proportion Explained   0.03807
## Cumulative Proportion       NA
## 
## Scaling 2 for species and site scores
## * Species are scaled proportional to eigenvalues
## * Sites are unscaled: weighted dispersion equal on all dimensions
## * General scaling constant of scores:  3.533199 
## 
## 
## Site scores (weighted sums of species scores)
## 
##      dbRDA1   dbRDA2   dbRDA3    dbRDA4   dbRDA5   dbRDA6
## 1  -0.07530 -0.62580 -0.63602 -2.643196 -0.34427 -23.8946
## 2   0.87242  0.77832 -0.32068  2.074749 -0.25800  41.3754
## 3  -0.81094  1.29200 -0.57735 -0.708420 -2.35800  -5.0035
## 4  -0.97450  0.91051  1.79849 -2.192492 -1.99248  -4.7134
## 5  -0.25889  0.74559 -1.30995 -1.418227  0.33185 -21.9720
## 6   0.41640  1.01888 -0.34551 -2.318324  1.61924 -24.8949
## 7  -0.98286 -0.67495  0.62017  1.174218  1.55155  -0.8974
## 8   0.93807  0.62635 -3.10878 -0.765073  2.41981  -3.0764
## 9  -0.62816 -0.77245 -1.30540  1.806170  1.30543 -41.2679
## 10 -1.03727 -0.50089  0.62572  0.369727  1.00378   4.7819
## 11  1.10743 -0.99448 -0.66515 -0.071830  0.83402 -14.8528
## 12  0.61624  0.76828  0.87110  0.357836  0.12125  38.2552
## 13 -0.96088 -0.72453  0.47977  1.569869  1.43399  -6.8692
## 14  1.41128  2.03453 -1.67256  0.776255  0.27648  29.9435
## 15 -1.00448 -0.59302  0.65010  0.833950  1.44643   3.4902
## 16  1.26382  0.81623  3.19230  1.380128  1.12523 -10.4959
## 17 -0.92726 -0.73637  0.33774  1.892276  1.07787  -8.2400
## 18 -1.04523 -1.01280 -0.59052  2.897871  0.35870 -14.2865
## 19  0.89987 -1.54255 -0.54435  1.738188 -2.61939  32.8867
## 20 -0.14608 -0.34079  0.78630  1.316880  1.70249  -9.1246
## 21 -1.30124 -1.46318 -1.27370  3.739798  0.91928 -16.7232
## 22 -0.07485  1.73636  0.30728 -1.760677 -0.46041   6.1504
## 23 -0.18008 -0.64105  0.55428 -1.310793 -0.82437  -3.6579
## 24  0.53054  0.08088  1.31623 -3.529128 -1.96164  41.6937
## 25 -0.50901 -0.82222 -1.26480  2.167458 -6.71178 -22.4535
## 26  0.54874  0.92050  1.12409  0.009258  0.07714  34.0773
## 27  1.25804 -0.69768  2.37255 -1.697654 -3.24936 -16.2349
## 28 -0.86554  0.32492 -0.68347  0.704689 -0.16329  13.0277
## 29  1.82276  0.98120  0.33891 -1.852448  4.33581  -8.6587
## 30  0.72974 -0.59907 -1.00943 -0.201449  1.77124  -4.5186
## 31 -1.02904  0.19515  1.05102 -1.468951 -1.36610  23.1359
## 32 -0.13877 -0.06050 -0.19569 -1.899083 -1.10403 -31.4856
## 33  1.01209 -0.59332  2.58400 -0.318337 -1.43566  39.0635
## 34 -0.39319 -1.79743  0.20497 -0.362972 -0.59693 -19.8004
## 35 -0.89506  1.46497 -0.05622 -1.876722 -1.86120   3.5844
## 36  0.87386  1.17335 -3.50567  1.196922  1.51036  -9.4970
## 37 -0.06267 -0.67495 -0.14975  0.389535  2.08495  11.1529
## 
## 
## Site constraints (linear combinations of constraining variables)
## 
##      dbRDA1   dbRDA2   dbRDA3   dbRDA4    dbRDA5   dbRDA6
## 1  -0.16605 -0.69736 -0.45052  0.39949 -1.140870  0.26766
## 2  -0.07683 -0.33692  0.61638  0.50606 -0.381728  1.66676
## 3  -0.10071  0.76299  0.51759  0.02910  0.003366 -0.26259
## 4  -1.55548  0.90013  0.66933 -0.51241 -0.415318 -0.13093
## 5  -0.41914  0.27994 -0.66104 -1.10909  0.384267  0.50060
## 6   0.63090  0.30135 -0.07724  0.37780  0.857007 -0.33310
## 7  -0.15252 -0.25647 -0.20458 -0.38652  0.792931  0.48033
## 8  -0.14874  0.12938 -0.39866 -0.78891  0.083814  0.27677
## 9  -0.14359  0.11628  0.95439 -0.04629  0.609992 -0.64671
## 10 -0.74154 -0.65088  0.87072 -1.17151 -0.875610 -1.73263
## 11  0.77571 -0.48314  0.02321 -0.76003  0.259283 -0.01494
## 12 -0.44387 -0.28442  0.54938 -0.01813  1.062736 -0.71676
## 13 -0.25608 -0.61494 -0.47847  0.74269  0.180828 -0.14346
## 14 -0.23636  0.67472 -0.87185  0.72759  0.046993 -0.61306
## 15 -0.75508  0.19781  0.41035  0.89449 -0.744271 -0.25040
## 16  0.97946  1.50039  0.74911  0.12808 -0.657115 -0.30363
## 17 -1.19679 -0.18091 -0.90174  0.16093 -0.229189  0.70454
## 18  0.23197  0.39183 -0.42932 -0.82808 -0.242760  0.64991
## 19  0.80868 -0.55920 -0.22435  0.32426 -0.445610 -0.51746
## 20 -0.19206  0.12003  0.24728  0.99541  0.747110  0.40576
## 21 -0.38477  0.09371 -0.03632  1.07969  1.383306 -0.25819
## 22  0.02352  0.34014  0.07462  0.03662 -0.352806  0.35626
## 23  0.18716 -0.91860 -0.02126 -0.69548  0.117329  0.60753
## 24  0.03851  0.57991 -0.70662 -1.10691  0.533700 -0.10095
## 25  0.61401 -0.01017 -0.51875  0.25795 -1.076732 -0.24901
## 26  0.42275  0.07490  0.99168  0.37950 -0.249480  0.19511
## 27  0.54478 -0.12356  0.31406  0.41569  0.237094  0.09772
## 28  0.02044  0.43106 -0.71781  0.33274 -0.428303 -0.26278
## 29  1.23235  0.15293  0.10463 -0.07285  0.107338  0.07885
## 30  0.02214 -1.33353 -0.37271  0.38762 -0.062960 -0.48751
## 31 -0.71637  0.11204 -0.85950  0.24734 -0.066602 -0.17455
## 32  0.01756 -0.71899 -0.27770  0.01727  0.106473 -0.73172
## 33  0.78296 -0.75806  0.27650 -0.27767 -0.031297 -0.01895
## 34 -0.33590 -0.63554  1.07989 -0.07860 -0.526645  0.97933
## 35 -0.15761  0.49842  0.66426 -0.05091  0.331527  0.71522
## 36  0.53360  0.88375 -0.84613  0.03571 -0.630019 -0.08456
## 37  0.31300  0.02099 -0.05877 -0.57261  0.712218  0.05155
## 
## 
## Biplot scores for constraining variables
## 
##               dbRDA1   dbRDA2  dbRDA3   dbRDA4  dbRDA5   dbRDA6
## T_av         0.63982  0.25527  0.1662 -0.34674 -0.2289 -0.48193
## O2_sat_av    0.02371 -0.02938 -0.8854 -0.18557  0.2035 -0.13356
## Con_av      -0.25315 -0.27391 -0.1407 -0.76998 -0.4264  0.07657
## COD_av      -0.13235  0.15034  0.4807  0.08407  0.1430 -0.57341
## NH4._av     -0.40051  0.14159  0.7569 -0.21438 -0.2051 -0.37005
## Nt_av       -0.31714 -0.23687  0.2952 -0.33830 -0.5347 -0.40917
## pool_riffle -0.34609 -0.15029 -0.4384  0.27287 -0.5547 -0.42271
## meander     -0.33563  0.52703 -0.3800 -0.21165 -0.5422 -0.13159
anova(spe.rda)
## Permutation test for dbrda under reduced model
## Permutation: free
## Number of permutations: 999
## 
## Model: dbrda(formula = spe.hel_bray ~ T_av + O2_sat_av + Con_av + COD_av + NH4._av + Nt_av + pool_riffle + meander, data = env_select)
##          Df SumOfSqs      F Pr(>F)
## Model     8   1.1664 1.2909  0.167
## Residual 28   3.1624
anova(spe.rda, by="term")
## Permutation test for dbrda under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## Model: dbrda(formula = spe.hel_bray ~ T_av + O2_sat_av + Con_av + COD_av + NH4._av + Nt_av + pool_riffle + meander, data = env_select)
##             Df SumOfSqs      F Pr(>F)  
## T_av         1   0.2734 2.4210  0.048 *
## O2_sat_av    1   0.1377 1.2195  0.302  
## Con_av       1   0.1613 1.4283  0.237  
## COD_av       1   0.0676 0.5990  0.695  
## NH4._av      1   0.1709 1.5131  0.212  
## Nt_av        1   0.0501 0.4439  0.787  
## pool_riffle  1   0.0657 0.5818  0.701  
## meander      1   0.2396 2.1210  0.070 .
## Residual    28   3.1624                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova.cca(spe.rda, step=1000);
## Permutation test for dbrda under reduced model
## Permutation: free
## Number of permutations: 999
## 
## Model: dbrda(formula = spe.hel_bray ~ T_av + O2_sat_av + Con_av + COD_av + NH4._av + Nt_av + pool_riffle + meander, data = env_select)
##          Df SumOfSqs      F Pr(>F)
## Model     8   1.1664 1.2909  0.164
## Residual 28   3.1624
anova.cca(spe.rda, step=1000, by="term");
## Permutation test for dbrda under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## Model: dbrda(formula = spe.hel_bray ~ T_av + O2_sat_av + Con_av + COD_av + NH4._av + Nt_av + pool_riffle + meander, data = env_select)
##             Df SumOfSqs      F Pr(>F)  
## T_av         1   0.2734 2.4210  0.057 .
## O2_sat_av    1   0.1377 1.2195  0.314  
## Con_av       1   0.1613 1.4283  0.245  
## COD_av       1   0.0676 0.5990  0.705  
## NH4._av      1   0.1709 1.5131  0.218  
## Nt_av        1   0.0501 0.4439  0.780  
## pool_riffle  1   0.0657 0.5818  0.699  
## meander      1   0.2396 2.1210  0.086 .
## Residual    28   3.1624                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
RsquareAdj(spe.rda)$adj.r.squared;
## [1] 0.06072755
RsquareAdj(spe.rda)$r.squared
## [1] 0.2694548

8.1.2 Effect of space on component community structure

# Same for spatial predictors
spe.rda <- dbrda(spe.hel_bray ~ netcen + updist, data = env_select)
anova(spe.rda)
## Permutation test for dbrda under reduced model
## Permutation: free
## Number of permutations: 999
## 
## Model: dbrda(formula = spe.hel_bray ~ netcen + updist, data = env_select)
##          Df SumOfSqs      F Pr(>F)  
## Model     2   0.5154 2.2975  0.036 *
## Residual 34   3.8135                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
RsquareAdj(spe.rda)$adj.r.squared
## [1] 0.06723433
mod0 <- dbrda(spe.hel_bray ~ 1, env_select[,c(9:10)])  # Model with intercept only  #edit_PH
mod1 <- dbrda(spe.hel_bray ~ ., env_select[,c(9:10)])  # Model with all explanatory variables  #edit_PH
step.res <- ordiR2step(mod0, mod1, direction = "both",perm.max = 200)
## Step: R2.adj= 0 
## Call: spe.hel_bray ~ 1 
##  
##                 R2.adjusted
## <All variables>  0.06723433
## + updist         0.04867127
## + netcen         0.04317133
## <none>           0.00000000
## 
##          Df    AIC      F Pr(>F)  
## + updist  1 54.314 2.8418  0.024 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: R2.adj= 0.04867127 
## Call: spe.hel_bray ~ updist 
##  
##                 R2.adjusted
## <All variables>  0.06723433
## + netcen         0.06723433
## <none>           0.04867127
## 
##          Df    AIC      F Pr(>F)
## + netcen  1 54.512 1.6965  0.152
step.res$anova  # Summary table
##                   R2.adj Df    AIC      F Pr(>F)  
## + updist        0.048671  1 54.314 2.8418  0.024 *
## <All variables> 0.067234                          
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
RsquareAdj(spe.rda)$adj.r.squared
## [1] 0.06723433
anova.cca(spe.rda, step=1000);
## Permutation test for dbrda under reduced model
## Permutation: free
## Number of permutations: 999
## 
## Model: dbrda(formula = spe.hel_bray ~ netcen + updist, data = env_select)
##          Df SumOfSqs      F Pr(>F)  
## Model     2   0.5154 2.2975  0.023 *
## Residual 34   3.8135                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova.cca(spe.rda, step=1000, by="term");
## Permutation test for dbrda under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## Model: dbrda(formula = spe.hel_bray ~ netcen + updist, data = env_select)
##          Df SumOfSqs      F Pr(>F)  
## netcen    1   0.3019 2.6920  0.033 *
## updist    1   0.2134 1.9029  0.113  
## Residual 34   3.8135                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
RsquareAdj(spe.rda)$adj.r.squared;
## [1] 0.06723433
RsquareAdj(spe.rda)$r.squared
## [1] 0.1190546

8.1.3 Variation partitioning

#Variation partitioning
spe.varpart1 <- varpart(spe.hel_bray, env_select[,1:8], env_select[,9:10])
plot(spe.varpart1,digits=2)

spe.varpart1
## 
## Partition of squared Bray distance in dbRDA 
## 
## Call: varpart(Y = spe.hel_bray, X = env_select[, 1:8], env_select[,
## 9:10])
## 
## Explanatory tables:
## X1:  env_select[, 1:8]
## X2:  env_select[, 9:10] 
## 
## No. of explanatory tables: 2 
## Total variation (SS): 4.3288 
## No. of observations: 37 
## 
## Partition table:
##                      Df R.squared Adj.R.squared Testable
## [a+c] = X1            8   0.26945       0.06073     TRUE
## [b+c] = X2            2   0.11905       0.06723     TRUE
## [a+b+c] = X1+X2      10   0.36716       0.12376     TRUE
## Individual fractions                                    
## [a] = X1|X2           8                 0.05653     TRUE
## [b] = X2|X1           2                 0.06303     TRUE
## [c]                   0                 0.00420    FALSE
## [d] = Residuals                         0.87624    FALSE
## ---
## Use function 'dbrda' to test significance of fractions of interest
anova.cca(dbrda(spe.hel_bray ~ T_av + O2_sat_av + Con_av + COD_av 
                     + NH4._av + Nt_av + pool_riffle + meander + Condition(netcen + updist),
                data=env_select), step=1000)
## Permutation test for dbrda under reduced model
## Permutation: free
## Number of permutations: 999
## 
## Model: dbrda(formula = spe.hel_bray ~ T_av + O2_sat_av + Con_av + COD_av + NH4._av + Nt_av + pool_riffle + meander + Condition(netcen + updist), data = env_select)
##          Df SumOfSqs      F Pr(>F)
## Model     8   1.0740 1.2742  0.195
## Residual 26   2.7395
anova.cca(dbrda(spe.hel_bray ~ netcen + updist+
                  Condition(T_av + O2_sat_av + Con_av + COD_av 
                     + NH4._av + Nt_av + pool_riffle + meander), data=env_select), step=1000)
## Permutation test for dbrda under reduced model
## Permutation: free
## Number of permutations: 999
## 
## Model: dbrda(formula = spe.hel_bray ~ netcen + updist + Condition(T_av + O2_sat_av + Con_av + COD_av + NH4._av + Nt_av + pool_riffle + meander), data = env_select)
##          Df SumOfSqs      F Pr(>F)  
## Model     2  0.42295 2.0071  0.047 *
## Residual 26  2.73945                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

8.2 Infra-communities

# Infracommunities: Bray-Curtis dissimilarities are calculated at the individual host level Hellinger-transformed parasite data and then averaged within site
# A dummy parasite species is added to avoid problems with non-infected fishes
data_infra <- na.omit(data[,c(1,22:24,26:32)])
data_infra_disp <- dispweight(data_infra[,-1])
braycurtis <- vegdist(decostand(cbind(data_infra_disp,rep(1,nrow(data_infra))), na.rm=T, method="hellinger"), method="bray", na.rm=T)
meandist_bray <- meandist(braycurtis, data_infra[,1])

# Check whether Euclidean and Bray-Curtis distances are comparable
braycurtis <- vegdist(decostand(cbind(data_infra_disp,rep(1,nrow(data_infra))), na.rm=T, method="hellinger"), method="bray", na.rm=T)
meandist_bray <- meandist(braycurtis, data_infra[,1])
euc <- vegdist(decostand(cbind(data_infra_disp,rep(1,nrow(data_infra))), na.rm=T, method="hellinger"), method="euc", na.rm=T)
meandist_euc <- meandist(euc, data_infra[,1])
plot(meandist_bray[1:37,1:37], meandist_euc[1:37,1:37])

mantel(meandist_bray[1:37,1:37], meandist_euc[1:37,1:37])
## 
## Mantel statistic based on Pearson's product-moment correlation 
## 
## Call:
## mantel(xdis = meandist_bray[1:37, 1:37], ydis = meandist_euc[1:37,      1:37]) 
## 
## Mantel statistic r: 0.9906 
##       Significance: 0.001 
## 
## Upper quantiles of permutations (null model):
##   90%   95% 97.5%   99% 
## 0.184 0.241 0.292 0.391 
## Permutation: free
## Number of permutations: 999
model_adonis = adonis(meandist_bray ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av 
                     + NH4._av + Nt_av + pool_riffle + meander + netcen + 
                       updist, data=environment2)
## 'adonis' will be deprecated: use 'adonis2' instead
model_adonis$aov.tab
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##             Df SumsOfSqs  MeanSqs F.Model      R2 Pr(>F)   
## avlength     1   0.05123 0.051226  6.1160 0.11570  0.006 **
## avcondition  1   0.00393 0.003927  0.4688 0.00887  0.622   
## T_av         1   0.00874 0.008741  1.0437 0.01974  0.369   
## O2_sat_av    1   0.02116 0.021158  2.5262 0.04779  0.105   
## Con_av       1   0.05861 0.058606  6.9971 0.13237  0.006 **
## COD_av       1   0.03543 0.035428  4.2298 0.08002  0.027 * 
## NH4._av      1   0.00275 0.002754  0.3288 0.00622  0.690   
## Nt_av        1   0.01049 0.010488  1.2521 0.02369  0.293   
## pool_riffle  1   0.00491 0.004914  0.5867 0.01110  0.563   
## meander      1   0.00959 0.009588  1.1447 0.02166  0.347   
## netcen       1   0.02319 0.023187  2.7684 0.05237  0.070 . 
## updist       1   0.01171 0.011707  1.3977 0.02644  0.283   
## Residuals   24   0.20102 0.008376         0.45403          
## Total       36   0.44274                  1.00000          
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# environmental variables
env_select <- environment2[,c("T_av", "O2_sat_av", "Con_av", "COD_av", "NH4._av", "Nt_av", "pool_riffle", "meander", "netcen", "updist")]
env_select$pool_riffle <- as.numeric(env_select$pool_riffle)
env_select$meander <- as.numeric(env_select$meander)

pca <- prcomp(env_select, scale.=T)
summary(pca)
## Importance of components:
##                           PC1    PC2    PC3    PC4     PC5     PC6     PC7
## Standard deviation     1.7124 1.5545 1.1221 1.0140 0.88807 0.79463 0.56647
## Proportion of Variance 0.2933 0.2416 0.1259 0.1028 0.07887 0.06314 0.03209
## Cumulative Proportion  0.2933 0.5349 0.6608 0.7636 0.84248 0.90563 0.93771
##                            PC8     PC9    PC10
## Standard deviation     0.50483 0.46939 0.38429
## Proportion of Variance 0.02549 0.02203 0.01477
## Cumulative Proportion  0.96320 0.98523 1.00000
plot(pca)

biplot(pca)

8.2.1 Effect of environment on infracommunity structure

# Assess the effect of environmental variables on parasite infracommunity dissimilarities using distance based RDA
spe.rda <- dbrda(meandist_bray ~ T_av + O2_sat_av + Con_av + COD_av 
                     + NH4._av + Nt_av + pool_riffle + meander, data = env_select)
anova(spe.rda)
## Permutation test for dbrda under reduced model
## Permutation: free
## Number of permutations: 999
## 
## Model: dbrda(formula = meandist_bray ~ T_av + O2_sat_av + Con_av + COD_av + NH4._av + Nt_av + pool_riffle + meander, data = env_select)
##          Df SumOfSqs      F Pr(>F)   
## Model     8  0.45036 1.3833  0.008 **
## Residual 28  1.13946                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
RsquareAdj(spe.rda)$adj.r.squared
## [1] 0.07850103
mod0 <- dbrda(meandist_bray ~ 1, env_select)  # Model with intercept only  #edit_PH
mod1 <- dbrda(meandist_bray ~ ., env_select)  # Model with all explanatory variables  #edit_PH
step.res <- ordiR2step(mod0, mod1, direction = "both",perm.max = 200)
## Step: R2.adj= 0 
## Call: meandist_bray ~ 1 
##  
##                   R2.adjusted
## <All variables>  0.0898822750
## + Con_av         0.0492715387
## + NH4._av        0.0376899022
## + Nt_av          0.0353578268
## + COD_av         0.0097867139
## + updist         0.0092890919
## + pool_riffle    0.0070398050
## + netcen         0.0034499960
## + O2_sat_av      0.0031240320
## <none>           0.0000000000
## + T_av          -0.0005031246
## + meander       -0.0037253892
## 
##          Df    AIC      F Pr(>F)   
## + Con_av  1 17.228 2.8657  0.004 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: R2.adj= 0.04927154 
## Call: meandist_bray ~ Con_av 
##  
##                 R2.adjusted
## <All variables>  0.08988227
## + COD_av         0.08426096
## + NH4._av        0.07058808
## + updist         0.06800179
## + O2_sat_av      0.06266812
## + netcen         0.05455352
## + meander        0.05428744
## + Nt_av          0.05273614
## <none>           0.04927154
## + pool_riffle    0.04923534
## + T_av           0.04694754
## 
##          Df    AIC      F Pr(>F)   
## + COD_av  1 16.768 2.3373  0.002 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: R2.adj= 0.08426096 
## Call: meandist_bray ~ Con_av + COD_av 
##  
##                 R2.adjusted
## + updist         0.09961037
## + meander        0.09248673
## <All variables>  0.08988227
## + netcen         0.08626581
## + pool_riffle    0.08559576
## <none>           0.08426096
## + T_av           0.08115910
## + O2_sat_av      0.08092067
## + Nt_av          0.07979657
## + NH4._av        0.07591794
step.res$anova  # Summary table
##                   R2.adj Df    AIC      F Pr(>F)   
## + Con_av        0.049272  1 17.228 2.8657  0.004 **
## + COD_av        0.084261  1 16.768 2.3373  0.002 **
## <All variables> 0.089882                           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
spe.rda <- dbrda(meandist_bray ~  Con_av + COD_av, env_select)
plot(spe.rda, scaling = 1) # it is for technical reasons not possible to plot both site and species scores

summary(spe.rda)
## 
## Call:
## dbrda(formula = meandist_bray ~ Con_av + COD_av, data = env_select) 
## 
## Partitioning of squared Unknown distance:
##               Inertia Proportion
## Total          1.5898     1.0000
## Constrained    0.2148     0.1351
## Unconstrained  1.3750     0.8649
## 
## Eigenvalues, and their contribution to the squared Unknown distance 
## 
## Importance of components:
##                       dbRDA1  dbRDA2   MDS1    MDS2   MDS3    MDS4    MDS5
## Eigenvalue            0.1863 0.02858 0.1972 0.15873 0.1207 0.08661 0.07076
## Proportion Explained  0.1172 0.01798 0.1240 0.09984 0.0759 0.05447 0.04451
## Cumulative Proportion     NA      NA     NA      NA     NA      NA      NA
##                          MDS6    MDS7    MDS8    MDS9   MDS10   MDS11   MDS12
## Eigenvalue            0.06231 0.04802 0.04613 0.04387 0.04231 0.03826 0.03703
## Proportion Explained  0.03919 0.03021 0.02901 0.02760 0.02661 0.02406 0.02329
## Cumulative Proportion      NA      NA      NA      NA      NA      NA      NA
##                         MDS13   MDS14   MDS15   MDS16   MDS17   MDS18   MDS19
## Eigenvalue            0.03557 0.03355 0.03200 0.02808 0.02797 0.02575 0.02487
## Proportion Explained  0.02237 0.02110 0.02013 0.01766 0.01759 0.01620 0.01564
## Cumulative Proportion      NA      NA      NA      NA      NA      NA      NA
##                         MDS20   MDS21   MDS22   MDS23   MDS24   MDS25   MDS26
## Eigenvalue            0.02377 0.02238 0.02182 0.02104 0.01990 0.01880 0.01782
## Proportion Explained  0.01495 0.01408 0.01372 0.01323 0.01252 0.01183 0.01121
## Cumulative Proportion      NA      NA      NA      NA      NA      NA      NA
##                         MDS27   MDS28    MDS29    MDS30    MDS31    MDS32
## Eigenvalue            0.01719 0.01622 0.011778 0.010499 0.007903 0.006958
## Proportion Explained  0.01081 0.01020 0.007408 0.006604 0.004971 0.004376
## Cumulative Proportion      NA      NA       NA       NA       NA       NA
##                          MDS33     iMDS1
## Eigenvalue            0.005670 -0.006426
## Proportion Explained  0.003566  0.004042
## Cumulative Proportion       NA        NA
## 
## Accumulated constrained eigenvalues
## Importance of components:
##                       dbRDA1  dbRDA2
## Eigenvalue            0.1863 0.02858
## Proportion Explained  0.8670 0.13304
## Cumulative Proportion 0.8670 1.00000
## 
## Scaling 2 for species and site scores
## * Species are scaled proportional to eigenvalues
## * Sites are unscaled: weighted dispersion equal on all dimensions
## * General scaling constant of scores:  2.750505 
## 
## 
## Site scores (weighted sums of species scores)
## 
##            dbRDA1   dbRDA2     MDS1     MDS2     MDS3      MDS4
## SITE 1  -0.215744  0.44381  0.01302  0.31819 -0.28116 -0.178171
## SITE 11 -0.507248  0.82144 -0.09719 -0.28198 -0.04679  0.342258
## SITE 12 -0.746029  0.15775  0.31672  0.25334 -0.08234 -0.151368
## SITE 13 -0.005847 -0.21986 -0.02027  0.91847 -0.12233 -0.332031
## SITE 14 -0.101987  0.28502  0.01314  0.55873 -0.09066 -0.223950
## SITE 15 -0.290401 -0.51714 -0.23614 -0.10408 -0.36670 -0.263022
## SITE 16  0.016467  0.22318  0.14693 -0.09879 -0.21010  0.156923
## SITE 17 -0.053594  0.78979  0.68552  0.22623  0.55763  0.068555
## SITE 18  0.741297 -0.34140  0.34156 -0.78781 -0.33039  0.372025
## SITE 19  2.204147  0.31903  0.07576 -0.77176 -0.32619  0.237586
## SITE 2  -0.334822  0.53901  0.20181  0.23760  0.25255  0.460075
## SITE 20  0.132501 -1.43824 -0.36307 -0.15683  0.91526  0.652586
## SITE 21 -0.518164  0.15765  0.51672  0.22690 -0.05180  0.071838
## SITE 22 -0.681351 -0.71637  0.44536  0.61813  0.24285 -0.440713
## SITE 23  0.085803 -0.07384 -0.18742 -0.58764 -0.50280  0.146049
## SITE 24 -0.506368 -1.41199 -0.02580  0.53865  0.25884  0.040947
## SITE 26  0.017028  0.61108  0.37526 -0.36715 -0.10492  0.425754
## SITE 28  0.121498  0.47609  0.59944 -0.22080 -0.07107  0.465237
## SITE 29 -0.014961  0.35365  0.04914 -0.10076 -0.08378  0.641604
## SITE 3  -0.711132 -0.09645  0.08469 -0.08718 -0.21517  0.038217
## SITE 30 -0.812113 -0.12175  0.55582 -0.16918 -0.25866  0.087617
## SITE 31  0.043616 -0.61571 -0.31459  0.47663 -0.16031 -0.651535
## SITE 32  0.562390  0.35603 -0.38737  0.06837 -0.20621 -0.072343
## SITE 33 -0.372173  0.03501  0.10319  0.85337 -0.01445 -0.328696
## SITE 34 -0.466651  0.48378  0.50738  0.02349 -0.08407  0.306511
## SITE 35  0.163917 -1.66724 -0.04520  0.22652  0.57157  0.198776
## SITE 36 -0.296043 -0.76361 -0.59812  0.40106 -0.04152 -0.249681
## SITE 38  0.302702  0.34509 -0.16274 -0.49724 -0.25040 -0.005416
## SITE 39 -0.978171  0.31142  0.36167  0.40699  0.02355  0.020399
## SITE 4  -0.141634  0.44343 -0.02774 -0.11080 -0.07834  0.180628
## SITE 40  0.496377 -0.01681 -0.38350 -0.04028 -0.50240 -0.167352
## SITE 41  1.117637  0.21808 -0.82799 -0.60061 -0.49878 -0.637014
## SITE 42  0.437517 -0.97421 -1.46550  0.08641  1.33505  0.880493
## SITE 5   0.522987  0.17736  0.27513  0.18757 -0.04206  0.152978
## SITE 6   0.423826 -0.07352 -0.91744 -0.10155 -0.42928 -0.789691
## SITE 7   0.013170  1.17661  0.57227 -1.15454  1.52741 -1.616432
## SITE 9   0.351552  0.32386 -0.18044 -0.38769 -0.23204  0.160360
## 
## 
## Site constraints (linear combinations of constraining variables)
## 
##            dbRDA1   dbRDA2     MDS1     MDS2     MDS3      MDS4
## SITE 1  -0.079483  0.44275  0.01302  0.31819 -0.28116 -0.178171
## SITE 11 -0.692558  0.69501 -0.09719 -0.28198 -0.04679  0.342258
## SITE 12 -0.387516 -0.29270  0.31672  0.25334 -0.08234 -0.151368
## SITE 13  0.692033  0.28074 -0.02027  0.91847 -0.12233 -0.332031
## SITE 14  0.273072  0.58517  0.01314  0.55873 -0.09066 -0.223950
## SITE 15 -0.503973 -0.49304 -0.23614 -0.10408 -0.36670 -0.263022
## SITE 16 -0.041630  0.31054  0.14693 -0.09879 -0.21010  0.156923
## SITE 17  0.286975  0.22462  0.68552  0.22623  0.55763  0.068555
## SITE 18  0.355899 -0.83587  0.34156 -0.78781 -0.33039  0.372025
## SITE 19  1.702755 -0.40556  0.07576 -0.77176 -0.32619  0.237586
## SITE 2  -0.068404  0.54101  0.20181  0.23760  0.25255  0.460075
## SITE 20  0.005618 -0.62143 -0.36307 -0.15683  0.91526  0.652586
## SITE 21 -0.085484 -0.08616  0.51672  0.22690 -0.05180  0.071838
## SITE 22 -0.015359 -0.63175  0.44536  0.61813  0.24285 -0.440713
## SITE 23 -0.475351 -0.06131 -0.18742 -0.58764 -0.50280  0.146049
## SITE 24 -0.138656 -0.63048 -0.02580  0.53865  0.25884  0.040947
## SITE 26 -0.086478  0.70999  0.37526 -0.36715 -0.10492  0.425754
## SITE 28  0.227529  0.47771  0.59944 -0.22080 -0.07107  0.465237
## SITE 29 -0.053389  0.16940  0.04914 -0.10076 -0.08378  0.641604
## SITE 3  -0.805011 -0.44279  0.08469 -0.08718 -0.21517  0.038217
## SITE 30 -0.711357 -0.97402  0.55582 -0.16918 -0.25866  0.087617
## SITE 31  0.197349 -0.41392 -0.31459  0.47663 -0.16031 -0.651535
## SITE 32  0.342870  0.50733 -0.38737  0.06837 -0.20621 -0.072343
## SITE 33  0.260419  0.25833  0.10319  0.85337 -0.01445 -0.328696
## SITE 34 -0.213627  0.34583  0.50738  0.02349 -0.08407  0.306511
## SITE 35  0.299656 -0.71657 -0.04520  0.22652  0.57157  0.198776
## SITE 36 -0.265626 -0.39128 -0.59812  0.40106 -0.04152 -0.249681
## SITE 38 -0.177490  0.11612 -0.16274 -0.49724 -0.25040 -0.005416
## SITE 39 -0.502932  0.21901  0.36167  0.40699  0.02355  0.020399
## SITE 4  -0.213933  0.31087 -0.02774 -0.11080 -0.07834  0.180628
## SITE 40  0.157625 -0.02847 -0.38350 -0.04028 -0.50240 -0.167352
## SITE 41  0.308647  0.03805 -0.82799 -0.60061 -0.49878 -0.637014
## SITE 42 -0.043511  0.13560 -1.46550  0.08641  1.33505  0.880493
## SITE 5   0.774214 -0.17770  0.27513  0.18757 -0.04206  0.152978
## SITE 6  -0.221478  0.23059 -0.91744 -0.10155 -0.42928 -0.789691
## SITE 7  -0.080948  0.35692  0.57227 -1.15454  1.52741 -1.616432
## SITE 9  -0.020466  0.24745 -0.18044 -0.38769 -0.23204  0.160360
## 
## 
## Biplot scores for constraining variables
## 
##        dbRDA1  dbRDA2 MDS1 MDS2 MDS3 MDS4
## Con_av 0.7628  0.6467    0    0    0    0
## COD_av 0.4413 -0.8974    0    0    0    0
anova(spe.rda)
## Permutation test for dbrda under reduced model
## Permutation: free
## Number of permutations: 999
## 
## Model: dbrda(formula = meandist_bray ~ Con_av + COD_av, data = env_select)
##          Df SumOfSqs      F Pr(>F)    
## Model     2  0.21484 2.6563  0.001 ***
## Residual 34  1.37498                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(spe.rda, by="term")
## Permutation test for dbrda under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## Model: dbrda(formula = meandist_bray ~ Con_av + COD_av, data = env_select)
##          Df SumOfSqs      F Pr(>F)    
## Con_av    1  0.12032 2.9752  0.001 ***
## COD_av    1  0.09452 2.3373  0.008 ** 
## Residual 34  1.37498                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

8.2.2 Effect of space on infracommunity structure

spe.rda <- dbrda(meandist_bray ~ netcen + updist, data = env_select)
anova(spe.rda)
## Permutation test for dbrda under reduced model
## Permutation: free
## Number of permutations: 999
## 
## Model: dbrda(formula = meandist_bray ~ netcen + updist, data = env_select)
##          Df SumOfSqs      F Pr(>F)
## Model     2  0.10742 1.2319  0.158
## Residual 34  1.48240
RsquareAdj(spe.rda)$adj.r.squared
## [1] 0.01271734
anova.cca(spe.rda, step=1000);
## Permutation test for dbrda under reduced model
## Permutation: free
## Number of permutations: 999
## 
## Model: dbrda(formula = meandist_bray ~ netcen + updist, data = env_select)
##          Df SumOfSqs      F Pr(>F)
## Model     2  0.10742 1.2319   0.15
## Residual 34  1.48240
anova.cca(spe.rda, step=1000, by="term");
## Permutation test for dbrda under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## Model: dbrda(formula = meandist_bray ~ netcen + updist, data = env_select)
##          Df SumOfSqs      F Pr(>F)
## netcen    1  0.04949 1.1352  0.264
## updist    1  0.05792 1.3285  0.135
## Residual 34  1.48240
RsquareAdj(spe.rda)$adj.r.squared;
## [1] 0.01271734
RsquareAdj(spe.rda)$r.squared
## [1] 0.06756637

8.2.3 Variation partitioning

#Variation partitioning
spe.varpart1 <- varpart(meandist_bray, env_select[,1:8], env_select[,9:10])
plot(spe.varpart1,digits=2)

spe.varpart1
## 
## Partition of squared Unknown user-supplied distance in dbRDA 
## 
## Call: varpart(Y = meandist_bray, X = env_select[, 1:8], env_select[,
## 9:10])
## 
## Explanatory tables:
## X1:  env_select[, 1:8]
## X2:  env_select[, 9:10] 
## 
## No. of explanatory tables: 2 
## Total variation (SS): 1.5898 
## No. of observations: 37 
## 
## Partition table:
##                      Df R.squared Adj.R.squared Testable
## [a+c] = X1            8   0.28328       0.07850     TRUE
## [b+c] = X2            2   0.06757       0.01272     TRUE
## [a+b+c] = X1+X2      10   0.34269       0.08988     TRUE
## Individual fractions                                    
## [a] = X1|X2           8                 0.07716     TRUE
## [b] = X2|X1           2                 0.01138     TRUE
## [c]                   0                 0.00134    FALSE
## [d] = Residuals                         0.91012    FALSE
## ---
## Use function 'dbrda' to test significance of fractions of interest
anova.cca(dbrda(meandist_bray ~ T_av + O2_sat_av + Con_av + COD_av 
                     + NH4._av + Nt_av + pool_riffle + meander + Condition(netcen + updist),
                data=env_select), step=1000)
## Permutation test for dbrda under reduced model
## Permutation: free
## Number of permutations: 999
## 
## Model: dbrda(formula = meandist_bray ~ T_av + O2_sat_av + Con_av + COD_av + NH4._av + Nt_av + pool_riffle + meander + Condition(netcen + updist), data = env_select)
##          Df SumOfSqs      F Pr(>F)  
## Model     8   0.4374 1.3603  0.014 *
## Residual 26   1.0450                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova.cca(dbrda(meandist_bray ~ netcen + updist+
                  Condition(T_av + O2_sat_av + Con_av + COD_av 
                     + NH4._av + Nt_av + pool_riffle + meander), data=env_select), step=1000)
## Permutation test for dbrda under reduced model
## Permutation: free
## Number of permutations: 999
## 
## Model: dbrda(formula = meandist_bray ~ netcen + updist + Condition(T_av + O2_sat_av + Con_av + COD_av + NH4._av + Nt_av + pool_riffle + meander), data = env_select)
##          Df SumOfSqs      F Pr(>F)
## Model     2  0.09446 1.1751  0.175
## Residual 26  1.04500